Line data Source code
1 : #ifdef ENABLE_DAGMC
2 :
3 : #include "MoabSkinner.h"
4 : #include "VariadicTable.h"
5 : #include "AuxiliarySystem.h"
6 : #include "BinUtility.h"
7 : #include "GeometryUtils.h"
8 : #include "UserErrorChecking.h"
9 : #include "DisplacedProblem.h"
10 : #include "MooseMeshElementConversionUtils.h"
11 :
12 : #include "libmesh/elem.h"
13 : #include "libmesh/enum_io_package.h"
14 : #include "libmesh/enum_order.h"
15 : #include "libmesh/enum_fe_family.h"
16 : #include "libmesh/equation_systems.h"
17 : #include "libmesh/system.h"
18 : #include "libmesh/mesh_tools.h"
19 :
20 : registerMooseObject("CardinalApp", MoabSkinner);
21 :
22 : InputParameters
23 148 : MoabSkinner::validParams()
24 : {
25 148 : InputParameters params = GeneralUserObject::validParams();
26 296 : params.addParam<bool>("verbose", false, "Whether to print diagnostic information");
27 :
28 : // temperature binning
29 296 : params.addRequiredParam<std::string>("temperature",
30 : "Temperature variable by which to bin elements");
31 444 : params.addRangeCheckedParam<Real>(
32 296 : "temperature_min", 0.0, "temperature_min >= 0.0", "Lower bound of temperature bins");
33 296 : params.addRequiredParam<Real>("temperature_max", "Upper bound of temperature bins");
34 296 : params.addRequiredRangeCheckedParam<unsigned int>(
35 : "n_temperature_bins", "n_temperature_bins > 0", "Number of temperature bins");
36 :
37 : // density binning
38 296 : params.addParam<std::string>("density", "Density variable by which to bin elements");
39 444 : params.addRangeCheckedParam<Real>(
40 296 : "density_min", 0.0, "density_min >= 0.0", "Lower bound of density bins");
41 296 : params.addParam<Real>("density_max", "Upper bound of density bins");
42 296 : params.addRangeCheckedParam<unsigned int>(
43 : "n_density_bins", "n_density_bins > 0", "Number of density bins");
44 296 : params.addParam<std::vector<SubdomainName>>(
45 : "material_blocks",
46 : "List of mesh subdomain names (or IDs) for which to assign material names in the generated "
47 : "DAGMC geometry. Must be provided together with 'material_names' where both have the same "
48 : "length and are listed in the same order."
49 : "Any subdomain not listed in 'material_blocks' will have no material assignment and OpenMC "
50 : "will"
51 : "default to void for that region. This parameter is optional when used in combination with "
52 : "OpenMCCellAverageProblem"
53 : "because in this case, the default is to auto-detect the materials in the original h5m file "
54 : "and re-apply them during skinning.");
55 296 : params.addParam<std::vector<std::string>>(
56 : "material_names",
57 : "Material names (or IDs) to assign to subdomains in the generated DAGMC geometry. "
58 : "Must be provided together with 'material_blocks' where both have the same length"
59 : "and listed in the same order. This parameter is optional when used in combination with "
60 : "OpenMCCellAverageProblem because in this case, the default is to auto-detect the materials "
61 : "in the original h5m file and re-apply them during skinning.");
62 444 : params.addRangeCheckedParam<Real>(
63 296 : "faceting_tol", 1e-4, "faceting_tol > 0", "Faceting tolerance for DagMC");
64 444 : params.addRangeCheckedParam<Real>(
65 296 : "geom_tol", 1e-6, "geom_tol > 0", "Geometry tolerance for DagMC");
66 :
67 296 : params.addParam<bool>(
68 296 : "build_graveyard", false, "Whether to build a graveyard around the geometry");
69 444 : params.addRangeCheckedParam<Real>(
70 : "graveyard_scale_inner",
71 296 : 1.01,
72 : "graveyard_scale_inner > 1",
73 : "Multiplier on mesh bounding box to form inner graveyard surface");
74 296 : params.addParam<Real>("graveyard_scale_outer",
75 296 : 1.10,
76 : "Multiplier on mesh bounding box to form outer graveyard surface");
77 296 : params.addParam<std::string>("implicit_complement_material",
78 : "Assigns OpenMC material name or ID to the implicit complement "
79 : "region. If not provided, void material is assigned by default.");
80 : // TODO: would be nice to support other file formats as well, like exodus
81 296 : params.addParam<bool>(
82 : "output_skins",
83 296 : false,
84 : "Whether the skinned MOAB mesh (skins generated from the "
85 : "libMesh [Mesh]) should be written to a file. The files will be named moab_skins_<n>.h5m, "
86 : "where <n> is the time step index. You can then visualize these files by running "
87 : "'mbconvert'.");
88 296 : params.addParam<bool>("output_full",
89 296 : false,
90 : "Whether the MOAB mesh (copied from the libMesh [Mesh]) should "
91 : "be written to a file. The files will be named moab_full_<n>.h5m, where "
92 : "<n> is the time step index. "
93 : "You can then visualize these files by running 'mbconvert'.");
94 296 : params.addParam<bool>("use_displaced_mesh",
95 296 : false,
96 : "Whether the skinned mesh should be generated from a displaced mesh ");
97 296 : params.addParam<std::vector<BoundaryName>>(
98 : "vacuum_bcs_surfaces",
99 : "Mesh sideset names or numeric sideset IDs to assign DAGMC vacuum boundary conditions to. "
100 : "Both string names and integer IDs are accepted. "
101 : "If not specified, the surface defaults to transmission.");
102 296 : params.addParam<std::vector<BoundaryName>>(
103 : "reflective_bcs_surfaces",
104 : "Mesh sideset names or numeric sideset IDs to assign DAGMC reflective boundary conditions "
105 : "to. Both string names and integer IDs are accepted. "
106 : "If not specified, the surface defaults to transmission.");
107 148 : params.addClassDescription("Re-generate the OpenMC geometry on-the-fly according to changes in "
108 : "the mesh geometry and/or contours in temperature and density");
109 148 : return params;
110 0 : }
111 :
112 74 : MoabSkinner::MoabSkinner(const InputParameters & parameters)
113 : : GeneralUserObject(parameters),
114 74 : _serialized_solution(NumericVector<Number>::build(_communicator).release()),
115 148 : _verbose(getParam<bool>("verbose")),
116 148 : _temperature_name(getParam<std::string>("temperature")),
117 148 : _temperature_min(getParam<Real>("temperature_min")),
118 148 : _temperature_max(getParam<Real>("temperature_max")),
119 148 : _n_temperature_bins(getParam<unsigned int>("n_temperature_bins")),
120 74 : _temperature_bin_width((_temperature_max - _temperature_min) / _n_temperature_bins),
121 148 : _bin_by_density(isParamValid("density")),
122 148 : _faceting_tol(getParam<Real>("faceting_tol")),
123 148 : _geom_tol(getParam<Real>("geom_tol")),
124 148 : _graveyard_scale_inner(getParam<double>("graveyard_scale_inner")),
125 148 : _graveyard_scale_outer(getParam<double>("graveyard_scale_outer")),
126 148 : _output_skins(getParam<bool>("output_skins")),
127 148 : _output_full(getParam<bool>("output_full")),
128 74 : _scaling(1.0),
129 74 : _n_write(0),
130 74 : _standalone(true),
131 74 : _tet_mesh_built(false),
132 342 : _set_bcs(isParamSetByUser("vacuum_bcs_surfaces") || isParamSetByUser("reflective_bcs_surfaces"))
133 : {
134 148 : _build_graveyard = getParam<bool>("build_graveyard");
135 148 : _use_displaced = getParam<bool>("use_displaced_mesh");
136 :
137 148 : if (isParamSetByUser("implicit_complement_material"))
138 : {
139 : // If the user specify a material that doesn't exist in materials.xml file, OpenMC
140 : // will catch the mistake.
141 2 : _set_implicit_complement_material = true;
142 : _implicit_complement_group_name =
143 8 : "mat:" + getParam<std::string>("implicit_complement_material") + "_comp";
144 : }
145 :
146 : // we can probably support this in the future, it's just not implemented yet
147 74 : if (!getMooseMesh().getMesh().is_serial())
148 0 : mooseError("MoabSkinner does not yet support distributed meshes!");
149 :
150 : // Create MOAB interface
151 74 : _moab = std::make_shared<moab::Core>();
152 :
153 : // Create a skinner
154 148 : skinner = std::make_unique<moab::Skinner>(_moab.get());
155 :
156 : // Create a geom topo tool
157 148 : gtt = std::make_unique<moab::GeomTopoTool>(_moab.get());
158 :
159 74 : if (_bin_by_density)
160 : {
161 48 : checkRequiredParam(parameters, "density_min", "binning by density");
162 48 : checkRequiredParam(parameters, "density_max", "binning by density");
163 48 : checkRequiredParam(parameters, "n_density_bins", "binning by density");
164 :
165 48 : _density_min = getParam<Real>("density_min");
166 48 : _density_max = getParam<Real>("density_max");
167 48 : _n_density_bins = getParam<unsigned int>("n_density_bins");
168 48 : _density_name = getParam<std::string>("density");
169 24 : _density_bin_width = (_density_max - _density_min) / _n_density_bins;
170 :
171 24 : if (_density_max < _density_min)
172 1 : paramError("density_max", "'density_max' must be greater than 'density_min'");
173 : }
174 : else
175 : {
176 100 : checkUnusedParam(parameters, "density_min", "not binning by density");
177 100 : checkUnusedParam(parameters, "density_max", "not binning by density");
178 100 : checkUnusedParam(parameters, "n_density_bins", "not binning by density");
179 :
180 50 : _n_density_bins = 1;
181 : }
182 :
183 73 : if (_build_graveyard)
184 : {
185 39 : if (_graveyard_scale_outer < _graveyard_scale_inner)
186 1 : paramError("graveyard_scale_outer",
187 : "'graveyard_scale_outer' must be greater than 'graveyard_scale_inner'!");
188 : }
189 : else
190 : {
191 68 : checkUnusedParam(parameters, "graveyard_scale_inner", "'build_graveyard' is false");
192 68 : checkUnusedParam(parameters, "graveyard_scale_outer", "'build_graveyard' is false");
193 : }
194 :
195 : // get variable numbers
196 72 : _temperature_var_num = getAuxiliaryVariableNumber(_temperature_name, "temperature");
197 70 : if (_bin_by_density)
198 : {
199 21 : if (_temperature_name == _density_name)
200 1 : mooseError("The 'temperature' and 'density' variables cannot be the same!");
201 :
202 20 : _density_var_num = getAuxiliaryVariableNumber(_density_name, "density");
203 : }
204 :
205 68 : if (_temperature_max <= _temperature_min)
206 1 : paramError("temperature_max", "'temperature_max' must be greater than 'temperature_min'");
207 :
208 305 : for (unsigned int i = 0; i < _n_temperature_bins + 1; ++i)
209 238 : _temperature_bin_bounds.push_back(_temperature_min + i * _temperature_bin_width);
210 :
211 263 : for (unsigned int i = 0; i < _n_density_bins + 1; ++i)
212 196 : _density_bin_bounds.push_back(_density_min + i * _density_bin_width);
213 :
214 : // node numberings for first-order tets
215 134 : _tet4_nodes.push_back({0, 1, 2, 3});
216 :
217 : // node numbers for second-order tets
218 134 : _tet10_nodes.push_back({0, 4, 6, 7});
219 134 : _tet10_nodes.push_back({1, 5, 4, 8});
220 134 : _tet10_nodes.push_back({2, 6, 5, 9});
221 134 : _tet10_nodes.push_back({7, 8, 9, 3});
222 134 : _tet10_nodes.push_back({4, 9, 7, 8});
223 134 : _tet10_nodes.push_back({4, 5, 9, 8});
224 134 : _tet10_nodes.push_back({4, 7, 9, 6});
225 134 : _tet10_nodes.push_back({4, 9, 5, 6});
226 :
227 67 : moab::MBErrorHandler_Init();
228 67 : }
229 :
230 : void
231 68 : MoabSkinner::finalize()
232 : {
233 68 : moab::MBErrorHandler_Finalize();
234 68 : }
235 :
236 : moab::ErrorCode
237 2123428 : MoabSkinner::check(const moab::ErrorCode input) const
238 : {
239 : #ifdef DEBUG
240 : MB_CHK_ERR(input);
241 : #endif
242 2123428 : return moab::MB_SUCCESS;
243 : }
244 :
245 : unsigned int
246 92 : MoabSkinner::getAuxiliaryVariableNumber(const std::string & name,
247 : const std::string & param_name) const
248 : {
249 92 : if (!_fe_problem.getAuxiliarySystem().hasVariable(name))
250 2 : paramError(param_name, "Cannot find auxiliary variable '", name, "'!");
251 :
252 : // we require these variables to be constant monomial
253 90 : auto type = _fe_problem.getAuxiliarySystem().getFieldVariable<Real>(0, name).feType();
254 90 : if (type.family != MONOMIAL || type.order != 0)
255 1 : paramError(param_name, "Auxiliary variable '", name, "' must be a CONSTANT MONOMIAL type!");
256 :
257 89 : return _fe_problem.getAuxiliarySystem().getFieldVariable<Real>(0, name).number();
258 : }
259 :
260 : MooseMesh &
261 121934 : MoabSkinner::getMooseMesh()
262 : {
263 156510 : if (_use_displaced && _fe_problem.getDisplacedProblem() == nullptr)
264 0 : mooseError("Displaced mesh was requested but the displaced problem does not exist. "
265 : "set use_displaced_mesh = False");
266 34576 : return ((_use_displaced && _fe_problem.getDisplacedProblem())
267 278444 : ? _fe_problem.getDisplacedProblem()->mesh()
268 209292 : : _fe_problem.mesh());
269 : }
270 :
271 : MeshBase &
272 83925 : MoabSkinner::getDAGMCGeometryMesh()
273 : {
274 83925 : if (_tet_mesh)
275 : return *_tet_mesh;
276 34519 : return getMooseMesh().getMesh();
277 : }
278 :
279 : void
280 136 : MoabSkinner::buildTetMesh()
281 : {
282 : _tet_mesh.reset();
283 136 : _tet_mesh_built = true;
284 :
285 136 : MeshBase & source = getMooseMesh().getMesh();
286 :
287 : // Decide whether conversion is needed and validate element types
288 : bool needs_conversion = false;
289 132268 : for (const auto * elem : source.active_element_ptr_range())
290 : {
291 66014 : if (elem->type() != TET4 && elem->type() != TET10)
292 : {
293 : needs_conversion = true;
294 : break;
295 : }
296 136 : }
297 :
298 136 : if (!needs_conversion)
299 120 : return;
300 :
301 16 : if (_verbose)
302 16 : _console << "MoabSkinner: non-tetrahedral elements detected. "
303 16 : << "Building internal TET4 copy of the mesh..." << std::endl;
304 :
305 32 : _tet_mesh = source.clone();
306 16 : _tet_mesh->all_first_order();
307 16 : _tet_mesh->prepare_for_use();
308 :
309 : // Stamp each element with its original ID; convert3DMeshToAllTet4 copies element extra
310 : // integers onto every child tet, so each tet can report the element it came from.
311 32 : const unsigned int orig_id_tag = _tet_mesh->add_elem_integer("orig_id");
312 : std::vector<std::pair<dof_id_type, bool>> elems_to_process;
313 13038 : for (auto * elem : _tet_mesh->active_element_ptr_range())
314 : {
315 13006 : elem->set_extra_integer(orig_id_tag, elem->id());
316 13006 : elems_to_process.emplace_back(elem->id(), true);
317 16 : }
318 :
319 : std::set<subdomain_id_type> sids;
320 16 : _tet_mesh->subdomain_ids(sids);
321 16 : const subdomain_id_type tmp_remove_sid = *sids.rbegin() + 1;
322 :
323 : // convert3DMeshToAllTet4() contracts and re-prepares the mesh internally.
324 : std::vector<dof_id_type> converted_ids;
325 16 : MooseMeshElementConversionUtils::convert3DMeshToAllTet4(
326 : *_tet_mesh, elems_to_process, converted_ids, tmp_remove_sid, true);
327 :
328 16 : if (_verbose)
329 16 : _console << "MoabSkinner: internal TET4 mesh has " << _tet_mesh->n_active_elem()
330 16 : << " elements (from " << source.n_active_elem() << " original elements)." << std::endl;
331 32 : }
332 :
333 : void
334 193 : MoabSkinner::initialize()
335 : {
336 193 : if (!_tet_mesh_built)
337 136 : buildTetMesh();
338 :
339 193 : findBlocks();
340 : _block_id_to_material_name.clear();
341 :
342 193 : const bool has_material_blocks = isParamValid("material_blocks");
343 193 : const bool has_material_names = isParamValid("material_names");
344 193 : if (has_material_blocks != has_material_names)
345 1 : mooseError("'material_blocks' and 'material_names' must both be provided if specifying "
346 : "materials for the skinner; you have provided only one of these parameters");
347 :
348 192 : if (has_material_blocks)
349 : {
350 57 : const auto & block_names = getParam<std::vector<SubdomainName>>("material_blocks");
351 114 : const auto & mat_names = getParam<std::vector<std::string>>("material_names");
352 57 : if (block_names.size() != mat_names.size())
353 4 : paramError("material_names",
354 2 : "'material_names' (length " + Moose::stringify(mat_names.size()) +
355 2 : ") must be the same length as 'material_blocks' (length " +
356 : Moose::stringify(block_names.size()) + ")");
357 :
358 : std::set<SubdomainID> seen_ids;
359 147 : for (const auto i : index_range(block_names))
360 : {
361 94 : if (!MooseMeshUtils::hasSubdomainName(getMooseMesh().getMesh(), block_names[i]))
362 2 : paramError("material_blocks",
363 1 : "Subdomain '" + std::string(block_names[i]) + "' not found in the mesh");
364 93 : auto id = getMooseMesh().getSubdomainID(block_names[i]);
365 : if (seen_ids.count(id))
366 : // A repeated entry for blocks (including if the same block is passed by its name then later
367 : // by its ID) should error here
368 2 : paramError("material_blocks",
369 1 : "Subdomain '" + std::string(block_names[i]) +
370 : "' was already listed in 'material_blocks' (either by name or ID). Each "
371 : "block can only be listed one time.");
372 92 : seen_ids.insert(id);
373 : auto it = _blocks.find(id);
374 92 : _block_id_to_material_name[it->first] = mat_names[i];
375 : }
376 : }
377 : else
378 : {
379 : // Neither 'material_blocks' nor 'material_names' was supplied in the input.
380 : // This is only valid when the skinner is controlled externally, e.g. an .h5m was
381 : // prepared outside of Cardinal and material names is set in setMaterialNames().
382 135 : if (_standalone)
383 1 : mooseError("When using the MoabSkinner independent of an OpenMCCellAverageProblem "
384 : "it requires both 'material_blocks' and 'material_names' parameters.");
385 600 : for (const auto & [subdomain_id, block_index] : _blocks)
386 466 : _block_id_to_material_name[subdomain_id] = _material_names[block_index];
387 : }
388 :
389 : std::vector<std::string> unassigned_blocks;
390 747 : for (const auto & [subdomain_id, block_index] : _blocks)
391 560 : if (_block_id_to_material_name.find(subdomain_id) == _block_id_to_material_name.end())
392 4 : unassigned_blocks.push_back(getMooseMesh().getSubdomainName(subdomain_id));
393 :
394 187 : if (!unassigned_blocks.empty())
395 4 : mooseWarning("The following mesh subdomains have no material assignment in MoabSkinner "
396 : "and will default to void in OpenMC:\n ",
397 7 : Moose::stringify(unassigned_blocks, "\n "));
398 :
399 : // Set spatial dimension in MOAB
400 186 : check(_moab->set_dimension(getMooseMesh().getMesh().spatial_dimension()));
401 :
402 : // Create a meshset representing all of the MOAB tets
403 186 : check(_moab->create_meshset(moab::MESHSET_SET, _all_tets));
404 :
405 186 : createTags();
406 :
407 186 : createMOABElems();
408 :
409 : // Resolve sideset names/IDs to boundary IDs and check that no boundary appears in both
410 : // 'vacuum_bcs_surfaces' and 'reflective_bcs_surfaces'
411 186 : if (_set_bcs)
412 : {
413 72 : if (isParamSetByUser("vacuum_bcs_surfaces"))
414 90 : _vacuum_bcs_surface_ids = boundaryNamesToIDs(
415 : getParam<std::vector<BoundaryName>>("vacuum_bcs_surfaces"), "vacuum_bcs_surfaces");
416 66 : if (isParamSetByUser("reflective_bcs_surfaces"))
417 : _reflective_bcs_surface_ids =
418 84 : boundaryNamesToIDs(getParam<std::vector<BoundaryName>>("reflective_bcs_surfaces"),
419 : "reflective_bcs_surfaces");
420 33 : checkBoundaryConditionOverlap();
421 : }
422 182 : }
423 :
424 : void
425 72 : MoabSkinner::execute()
426 : {
427 72 : if (_standalone)
428 23 : update();
429 68 : }
430 :
431 : void
432 39 : MoabSkinner::setUseDisplacedMesh(const bool & use)
433 : {
434 39 : if ((use != _use_displaced) && isParamSetByUser("use_displaced_mesh"))
435 0 : mooseWarning("Overriding 'use_displaced_mesh' to " + std::to_string(use) +
436 : " to match the displaced problem action.");
437 39 : _use_displaced = use;
438 39 : }
439 :
440 : void
441 76 : MoabSkinner::update()
442 : {
443 76 : _console << "Skinning geometry into " << _n_temperature_bins << " temperature bins, "
444 76 : << _n_density_bins << " density bins, and " << _n_block_bins << " block bins... "
445 76 : << std::endl;
446 :
447 76 : if (_use_displaced && _standalone)
448 : {
449 : // we are responsible for updating the mesh if running in standalone mode; otherwise, the
450 : // OpenMCCellAverageProblem class does it
451 4 : _fe_problem.getDisplacedProblem()->updateMesh();
452 : }
453 :
454 : // Clear MOAB mesh data from last timestep
455 76 : reset();
456 :
457 76 : _serialized_solution->init(_fe_problem.getAuxiliarySystem().sys().n_dofs(), false, SERIAL);
458 76 : _fe_problem.getAuxiliarySystem().solution().localize(*_serialized_solution);
459 :
460 : // Re-initialise the mesh data
461 76 : initialize();
462 :
463 202 : if (isParamValid("material_blocks") && isParamValid("material_names"))
464 25 : _console << "MoabSkinner updating material assignments..." << std::endl;
465 :
466 : // Sort libMesh elements into bins
467 76 : sortElemsByResults();
468 :
469 : // Find the surfaces of local temperature regions
470 72 : findSurfaces();
471 68 : }
472 :
473 : void
474 193 : MoabSkinner::findBlocks()
475 : {
476 : _blocks.clear();
477 :
478 : int i = 0;
479 764 : for (const auto & b : getMooseMesh().meshSubdomains())
480 571 : _blocks[b] = i++;
481 :
482 193 : _n_block_bins = _blocks.size();
483 193 : }
484 :
485 : void
486 186 : MoabSkinner::createMOABElems()
487 : {
488 : // Clear prior results
489 : _id_to_elem_handles.clear();
490 : _node_id_to_handle.clear();
491 : _elem_handle_to_id.clear();
492 :
493 : double coords[3];
494 :
495 186 : MeshBase & geom_mesh = getDAGMCGeometryMesh();
496 :
497 : // Save all the node information
498 55404 : for (const auto & node : geom_mesh.node_ptr_range())
499 : {
500 : // Fetch coords (and scale to correct units)
501 55032 : coords[0] = _scaling * (*node)(0);
502 55032 : coords[1] = _scaling * (*node)(1);
503 55032 : coords[2] = _scaling * (*node)(2);
504 :
505 : // Add node to MOAB database and get handle
506 55032 : moab::EntityHandle ent(0);
507 55032 : check(_moab->create_vertex(coords, ent));
508 :
509 : // Save mapping of libMesh IDs to MOAB vertex handles
510 55032 : _node_id_to_handle[node->id()] = ent;
511 186 : }
512 :
513 : moab::Range all_elems;
514 :
515 : // Iterate over elements in the mesh
516 380530 : for (const auto & elem : geom_mesh.active_element_ptr_range())
517 : {
518 190079 : auto nodeSets = getTetSets(elem->type());
519 :
520 : // Get the connectivity
521 : std::vector<dof_id_type> conn_libmesh;
522 190079 : elem->connectivity(0, libMesh::IOPackage::VTK, conn_libmesh);
523 :
524 : // Loop over sub tets
525 402782 : for (const auto & nodeSet : nodeSets)
526 : {
527 : // Set MOAB connectivity
528 212703 : std::vector<moab::EntityHandle> conn(NODES_PER_MOAB_TET);
529 1063515 : for (unsigned int i = 0; i < NODES_PER_MOAB_TET; ++i)
530 : {
531 : // Get the elem node index of the ith node of the sub-tet
532 850812 : unsigned int nodeIndex = nodeSet.at(i);
533 1701624 : conn[i] = _node_id_to_handle[conn_libmesh.at(nodeIndex)];
534 : }
535 :
536 : // Create an element in MOAB database
537 212703 : moab::EntityHandle ent(0);
538 212703 : check(_moab->create_element(moab::MBTET, conn.data(), NODES_PER_MOAB_TET, ent));
539 :
540 : // Save mapping between libMesh ids and moab handles
541 212703 : auto id = elem->id();
542 212703 : if (_id_to_elem_handles.find(id) == _id_to_elem_handles.end())
543 380158 : _id_to_elem_handles[id] = std::vector<moab::EntityHandle>();
544 :
545 212703 : _id_to_elem_handles[id].push_back(ent);
546 212703 : _elem_handle_to_id[ent] = id;
547 :
548 : // Save the handle for adding to entity sets
549 212703 : all_elems.insert(ent);
550 212703 : }
551 190265 : }
552 :
553 : // Add the elems to the full meshset
554 186 : check(_moab->add_entities(_all_tets, all_elems));
555 :
556 : // Save the first elem
557 186 : offset = all_elems.front();
558 186 : }
559 :
560 : const std::vector<std::vector<unsigned int>> &
561 190079 : MoabSkinner::getTetSets(ElemType type) const
562 : {
563 190079 : if (type == TET4)
564 186847 : return _tet4_nodes;
565 3232 : else if (type == TET10)
566 3232 : return _tet10_nodes;
567 : else
568 0 : mooseError("The MoabSkinner can only be used with a tetrahedral [Mesh]! If your mesh "
569 : "contains other element types, MoabSkinner should convert them automatically.");
570 : }
571 :
572 : void
573 186 : MoabSkinner::createTags()
574 : {
575 : // First some built-in MOAB tag types
576 186 : check(_moab->tag_get_handle(GEOM_DIMENSION_TAG_NAME,
577 : 1,
578 : moab::MB_TYPE_INTEGER,
579 186 : geometry_dimension_tag,
580 : moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
581 :
582 186 : check(_moab->tag_get_handle(GLOBAL_ID_TAG_NAME,
583 : 1,
584 : moab::MB_TYPE_INTEGER,
585 186 : id_tag,
586 : moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
587 :
588 186 : check(_moab->tag_get_handle(CATEGORY_TAG_NAME,
589 : CATEGORY_TAG_SIZE,
590 : moab::MB_TYPE_OPAQUE,
591 186 : category_tag,
592 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
593 :
594 186 : check(_moab->tag_get_handle(NAME_TAG_NAME,
595 : NAME_TAG_SIZE,
596 : moab::MB_TYPE_OPAQUE,
597 186 : name_tag,
598 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
599 :
600 : // Some tags needed for DagMC
601 186 : check(_moab->tag_get_handle("FACETING_TOL",
602 : 1,
603 : moab::MB_TYPE_DOUBLE,
604 186 : faceting_tol_tag,
605 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
606 :
607 186 : check(_moab->tag_get_handle("GEOMETRY_RESABS",
608 : 1,
609 : moab::MB_TYPE_DOUBLE,
610 186 : geometry_resabs_tag,
611 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
612 :
613 : // Set the values for DagMC faceting / geometry tolerance tags on the mesh entity set
614 186 : check(_moab->tag_set_data(faceting_tol_tag, &_all_tets, 1, &_faceting_tol));
615 186 : check(_moab->tag_set_data(geometry_resabs_tag, &_all_tets, 1, &_geom_tol));
616 186 : }
617 :
618 : void
619 1121 : MoabSkinner::createGroup(const unsigned int & id,
620 : const std::string & name,
621 : moab::EntityHandle & group_set)
622 : {
623 1121 : check(_moab->create_meshset(moab::MESHSET_SET, group_set));
624 2242 : setTags(group_set, name, "Group", id, 4);
625 1121 : }
626 :
627 : void
628 1895 : MoabSkinner::createVol(const unsigned int & id,
629 : moab::EntityHandle & volume_set,
630 : moab::EntityHandle group_set)
631 : {
632 1895 : check(_moab->create_meshset(moab::MESHSET_SET, volume_set));
633 :
634 3790 : setTags(volume_set, "", "Volume", id, 3);
635 :
636 : // Add the volume to group
637 1895 : check(_moab->add_entities(group_set, &volume_set, 1));
638 1895 : }
639 :
640 : void
641 4371 : MoabSkinner::createSurf(const unsigned int & id,
642 : moab::EntityHandle & surface_set,
643 : moab::Range & faces,
644 : const std::vector<VolData> & voldata,
645 : BoundaryConditionType bc_type)
646 : {
647 : // Create meshset
648 4371 : check(_moab->create_meshset(moab::MESHSET_SET, surface_set));
649 :
650 : // Set tags
651 8742 : setTags(surface_set, "", "Surface", id, 2);
652 :
653 : // Add tris to the surface
654 4371 : check(_moab->add_entities(surface_set, faces));
655 :
656 : // Create entry in map
657 8742 : surfsToVols[surface_set] = std::vector<VolData>();
658 :
659 : // Add volume to list associated with this surface
660 11156 : for (const auto & data : voldata)
661 6785 : updateSurfData(surface_set, data);
662 :
663 4371 : recordBoundaryConditionSurface(surface_set, bc_type);
664 4371 : }
665 :
666 : void
667 7649 : MoabSkinner::updateSurfData(moab::EntityHandle surface_set, const VolData & data)
668 : {
669 : // Add the surface to the volume set
670 7649 : check(_moab->add_parent_child(data.vol, surface_set));
671 :
672 : // Set the surfaces sense
673 7649 : gtt->set_sense(surface_set, data.vol, int(data.sense));
674 :
675 7649 : surfsToVols[surface_set].push_back(data);
676 7649 : }
677 :
678 : void
679 7387 : MoabSkinner::setTags(
680 : moab::EntityHandle ent, std::string name, std::string category, unsigned int id, int dim)
681 : {
682 : // Set the name tag
683 7387 : if (name != "")
684 2242 : setTagData(name_tag, ent, name, NAME_TAG_SIZE);
685 :
686 : // Set the category tag
687 7387 : if (category != "")
688 14774 : setTagData(category_tag, ent, category, CATEGORY_TAG_SIZE);
689 :
690 : // Set the dimension tag
691 7387 : setTagData(geometry_dimension_tag, ent, &dim);
692 :
693 : // Set the id tag
694 7387 : setTagData(id_tag, ent, &id);
695 7387 : }
696 :
697 : void
698 8508 : MoabSkinner::setTagData(moab::Tag tag, moab::EntityHandle ent, std::string data, unsigned int SIZE)
699 : {
700 8508 : auto namebuf = new char[SIZE];
701 8508 : memset(namebuf, '\0', SIZE); // fill C char array with null
702 8508 : strncpy(namebuf, data.c_str(), SIZE - 1);
703 8508 : check(_moab->tag_set_data(tag, &ent, 1, namebuf));
704 8508 : delete[] namebuf;
705 8508 : }
706 :
707 : void
708 14774 : MoabSkinner::setTagData(moab::Tag tag, moab::EntityHandle ent, void * data)
709 : {
710 14774 : check(_moab->tag_set_data(tag, &ent, 1, data));
711 14774 : }
712 :
713 : unsigned int
714 134 : MoabSkinner::nBins() const
715 : {
716 134 : return _n_block_bins * _n_density_bins * _n_temperature_bins;
717 : }
718 :
719 : void
720 76 : MoabSkinner::sortElemsByResults()
721 : {
722 76 : _elem_bins.clear();
723 76 : _elem_bins.resize(nBins());
724 :
725 : // accumulate information for printing diagnostics
726 76 : std::vector<unsigned int> n_block_hits(_n_block_bins, 0);
727 76 : std::vector<unsigned int> n_temp_hits(_n_temperature_bins, 0);
728 76 : std::vector<unsigned int> n_density_hits(_n_density_bins, 0);
729 :
730 76 : if (!_tet_mesh)
731 : {
732 34523 : for (unsigned int e = 0; e < getMooseMesh().nElem(); ++e)
733 : {
734 34461 : const Elem * const elem = getMooseMesh().queryElemPtr(e);
735 34461 : if (!elem)
736 0 : continue;
737 :
738 34461 : auto iMat = getSubdomainBin(elem);
739 34461 : n_block_hits[iMat] += 1;
740 :
741 34461 : auto iDenBin = getDensityBin(elem);
742 34459 : n_density_hits[iDenBin] += 1;
743 :
744 34459 : auto iBin = getTemperatureBin(elem);
745 34457 : n_temp_hits[iBin] += 1;
746 :
747 34457 : _elem_bins.at(getBin(iBin, iDenBin, iMat)).insert(elem->id());
748 : }
749 : }
750 : else
751 : {
752 : // aux vars live on the original mesh; _elem_bins
753 : // must hold tet-mesh IDs (keyed the same as _id_to_elem_handles).
754 :
755 : // compute bin index for every original element.
756 : std::unordered_map<dof_id_type, unsigned int> orig_id_to_bin;
757 :
758 8638 : for (unsigned int e = 0; e < getMooseMesh().nElem(); ++e)
759 : {
760 8628 : const Elem * const elem = getMooseMesh().queryElemPtr(e);
761 8628 : if (!elem)
762 0 : continue;
763 :
764 8628 : auto iMat = getSubdomainBin(elem);
765 8628 : n_block_hits[iMat] += 1;
766 :
767 8628 : auto iDenBin = getDensityBin(elem);
768 8628 : n_density_hits[iDenBin] += 1;
769 :
770 8628 : auto iBin = getTemperatureBin(elem);
771 8628 : n_temp_hits[iBin] += 1;
772 :
773 8628 : orig_id_to_bin[elem->id()] = getBin(iBin, iDenBin, iMat);
774 : }
775 :
776 : // insert each tet ID into its parent's bin.
777 10 : const unsigned int orig_id_tag = _tet_mesh->get_elem_integer_index("orig_id");
778 49388 : for (const auto * tet : _tet_mesh->active_element_ptr_range())
779 : {
780 49368 : const dof_id_type orig_id = tet->get_extra_integer(orig_id_tag);
781 :
782 : auto bin_it = orig_id_to_bin.find(orig_id);
783 49368 : if (bin_it == orig_id_to_bin.end())
784 0 : mooseError("MoabSkinner::sortElemsByResults(): original elem id=",
785 : orig_id,
786 : " not found in bin map");
787 :
788 49368 : _elem_bins.at(bin_it->second).insert(tet->id());
789 10 : }
790 : }
791 :
792 72 : if (_verbose)
793 : {
794 69 : VariadicTable<unsigned int, std::string, unsigned int> vtt({"Bin", "Range (K)", "# Elems"});
795 69 : VariadicTable<unsigned int, std::string, unsigned int> vtd({"Bin", "Range (kg/m3)", "# Elems"});
796 :
797 266 : for (unsigned int i = 0; i < _n_temperature_bins; ++i)
798 197 : vtt.addRow(i,
799 591 : std::to_string(_temperature_bin_bounds[i]) + " to " +
800 394 : std::to_string(_temperature_bin_bounds[i + 1]),
801 197 : n_temp_hits[i]);
802 :
803 232 : for (unsigned int i = 0; i < _n_density_bins; ++i)
804 163 : vtd.addRow(i,
805 489 : std::to_string(_density_bin_bounds[i]) + " to " +
806 326 : std::to_string(_density_bin_bounds[i + 1]),
807 163 : n_density_hits[i]);
808 :
809 69 : _console << "\nMapping of Elements to Temperature Bins:" << std::endl;
810 69 : vtt.print(_console);
811 69 : _console << std::endl;
812 :
813 69 : if (_bin_by_density)
814 : {
815 21 : _console << "\n\nMapping of Elements to Density Bins:" << std::endl;
816 21 : vtd.print(_console);
817 21 : _console << std::endl;
818 : }
819 69 : }
820 72 : }
821 :
822 : unsigned int
823 100468 : MoabSkinner::getTemperatureBin(const Elem * const elem) const
824 : {
825 100468 : auto dof = elem->dof_number(_fe_problem.getAuxiliarySystem().number(), _temperature_var_num, 0);
826 100468 : auto value = (*_serialized_solution)(dof);
827 :
828 : // TODO: add option to truncate instead
829 100468 : if ((_temperature_min - value) > BIN_TOLERANCE)
830 1 : mooseError("Variable '",
831 1 : _temperature_name,
832 : "' has value below minimum range of bins. "
833 : "Please decrease 'temperature_min'.\n\n"
834 : " value: ",
835 : value,
836 : "\n temperature_min: ",
837 : _temperature_min);
838 :
839 100467 : if ((value - _temperature_max) > BIN_TOLERANCE)
840 1 : mooseError("Variable '",
841 1 : _temperature_name,
842 : "' has value above maximum range of bins. "
843 : "Please increase 'temperature_max'.\n\n"
844 : " value: ",
845 : value,
846 : "\n temperature_max: ",
847 : _temperature_max);
848 :
849 100466 : return bin_utility::linearBin(value, _temperature_bin_bounds);
850 : }
851 :
852 : unsigned int
853 98604 : MoabSkinner::getDensityBin(const Elem * const elem) const
854 : {
855 98604 : if (!_bin_by_density)
856 : return 0;
857 :
858 56106 : auto dof = elem->dof_number(_fe_problem.getAuxiliarySystem().number(), _density_var_num, 0);
859 56106 : auto value = (*_serialized_solution)(dof);
860 :
861 : // TODO: add option to truncate instead
862 56106 : if ((_density_min - value) > BIN_TOLERANCE)
863 1 : mooseError("Variable '",
864 1 : _density_name,
865 : "' has value below minimum range of bins. "
866 : "Please decrease 'density_min'.\n\n"
867 : " value: ",
868 : value,
869 : "\n density_min: ",
870 1 : _density_min);
871 :
872 56105 : if ((value - _density_max) > BIN_TOLERANCE)
873 1 : mooseError("Variable '",
874 1 : _density_name,
875 : "' has value above maximum range of bins. "
876 : "Please increase 'density_max'.\n\n"
877 : " value: ",
878 : value,
879 : "\n density_max: ",
880 1 : _density_max);
881 :
882 56104 : return bin_utility::linearBin(value, _density_bin_bounds);
883 : }
884 :
885 : std::string
886 1057 : MoabSkinner::materialName(const unsigned int & block,
887 : const unsigned int & density,
888 : const unsigned int & temp) const
889 : {
890 2385 : for (const auto & [subdomain_id, block_index] : _blocks)
891 : {
892 2385 : if (block_index == block)
893 : {
894 : auto it = _block_id_to_material_name.find(subdomain_id);
895 1057 : if (it != _block_id_to_material_name.end())
896 1053 : return "mat:" + it->second;
897 :
898 4 : return "mat:void";
899 : }
900 : }
901 : mooseAssert(false, "could not find subdomain for block index " + Moose::stringify(block));
902 0 : return "";
903 : }
904 :
905 : std::set<BoundaryID>
906 59 : MoabSkinner::boundaryNamesToIDs(const std::vector<BoundaryName> & names,
907 : const std::string & param_name)
908 : {
909 : std::set<BoundaryID> ids;
910 223 : for (const auto & name : names)
911 : {
912 167 : const BoundaryID id = getMooseMesh().getBoundaryID(name);
913 167 : if (!getMooseMesh().meshSidesetIds().count(id))
914 : {
915 3 : if (getMooseMesh().meshNodesetIds().count(id))
916 1 : paramError(param_name,
917 : "Boundary '",
918 : name,
919 : "' is a nodeset; boundary conditions can only be assigned to sidesets.");
920 : else
921 2 : paramError(param_name, "Boundary '", name, "' does not exist in the mesh.");
922 : }
923 164 : ids.insert(id);
924 : }
925 56 : return ids;
926 : }
927 :
928 : void
929 33 : MoabSkinner::checkBoundaryConditionOverlap() const
930 : {
931 : std::set<BoundaryID> overlap;
932 33 : std::set_intersection(_vacuum_bcs_surface_ids.begin(),
933 : _vacuum_bcs_surface_ids.end(),
934 : _reflective_bcs_surface_ids.begin(),
935 : _reflective_bcs_surface_ids.end(),
936 : std::inserter(overlap, overlap.begin()));
937 33 : if (overlap.empty())
938 32 : return;
939 :
940 : std::string ids_str;
941 2 : for (const auto id : overlap)
942 2 : ids_str += (ids_str.empty() ? "" : ", ") + std::to_string(id);
943 :
944 1 : paramError("vacuum_bcs_surfaces",
945 : "The following sideset ID(s) appear in both 'vacuum_bcs_surfaces' and "
946 : "'reflective_bcs_surfaces': ",
947 : ids_str);
948 : }
949 :
950 : MoabSkinner::BoundaryConditionType
951 25357 : MoabSkinner::boundaryConditionType(const Elem * const elem,
952 : const unsigned int side,
953 : const libMesh::BoundaryInfo & boundary_info) const
954 : {
955 : std::vector<boundary_id_type> side_bids;
956 25357 : boundary_info.boundary_ids(elem, side, side_bids);
957 :
958 : std::vector<boundary_id_type> vacuum_bids, reflective_bids;
959 31082 : for (const auto bid : side_bids)
960 : {
961 : if (_vacuum_bcs_surface_ids.count(bid))
962 657 : vacuum_bids.push_back(bid);
963 : if (_reflective_bcs_surface_ids.count(bid))
964 2321 : reflective_bids.push_back(bid);
965 : }
966 :
967 25357 : if (!vacuum_bids.empty() && !reflective_bids.empty())
968 : {
969 : // Assemble sideset names (resolves to IDs) for a helpful error message
970 2 : auto names = [&boundary_info](const std::vector<boundary_id_type> & bids)
971 : {
972 : std::string out;
973 4 : for (const auto bid : bids)
974 : {
975 2 : const auto & name = boundary_info.get_sideset_name(bid);
976 2 : out += (out.empty() ? "" : ", ");
977 4 : out += "'" + (name.empty() ? std::to_string(bid) : name) + "'";
978 : }
979 2 : return out;
980 1 : };
981 :
982 2 : mooseError("Element ",
983 1 : elem->id(),
984 : ", side ",
985 : side,
986 : " is assigned both vacuum and reflective boundary conditions, because this side "
987 : "belongs to sideset(s) ",
988 1 : names(vacuum_bids),
989 : " (vacuum) and ",
990 1 : names(reflective_bids),
991 : " (reflective). An element side cannot be assigned competing types of surface "
992 : "boundary conditions.");
993 : }
994 :
995 25356 : if (!vacuum_bids.empty())
996 : return BoundaryConditionType::Vacuum;
997 24700 : if (!reflective_bids.empty())
998 2320 : return BoundaryConditionType::Reflective;
999 : return BoundaryConditionType::Transmission;
1000 25356 : }
1001 :
1002 : void
1003 3693 : MoabSkinner::splitSkinByBoundaryCondition(const moab::Range & region,
1004 : const moab::Range & skin,
1005 : moab::Range & transmission_tris,
1006 : moab::Range & vacuum_tris,
1007 : moab::Range & reflective_tris)
1008 : {
1009 : // all triangles have transmission BC until sorted
1010 3693 : transmission_tris = skin;
1011 :
1012 3693 : if (!_set_bcs || skin.empty())
1013 : return;
1014 :
1015 34 : MeshBase & geom_mesh = getDAGMCGeometryMesh();
1016 : const auto & boundary_info = geom_mesh.get_boundary_info();
1017 :
1018 : // Number of entities passed to get_adjacencies per call
1019 : constexpr int num_entities = 1;
1020 : // Topological dimension of the faces we want (2 = triangles for a 3D tet mesh)
1021 : constexpr int surface_dimension = 2;
1022 :
1023 12713 : for (const auto tet : region)
1024 : {
1025 : const auto elem_id_it = _elem_handle_to_id.find(tet);
1026 : mooseAssert(elem_id_it != _elem_handle_to_id.end(),
1027 : "Could not map a MOAB tet back to a libMesh element while assigning "
1028 : "DAGMC boundary conditions.");
1029 :
1030 6340 : const Elem * const elem = geom_mesh.query_elem_ptr(elem_id_it->second);
1031 : mooseAssert(elem,
1032 : "Could not find libMesh element " + std::to_string(elem_id_it->second) +
1033 : " while assigning DAGMC boundary conditions.");
1034 :
1035 31696 : for (const auto side : make_range(elem->n_sides()))
1036 : {
1037 25357 : const auto bc_type = boundaryConditionType(elem, side, boundary_info);
1038 25356 : if (bc_type == BoundaryConditionType::Transmission)
1039 22380 : continue;
1040 :
1041 : // Build the set of MOAB vertex handles covering this libMesh side
1042 2976 : std::unique_ptr<const Elem> side_elem = elem->build_side_ptr(side);
1043 : std::set<moab::EntityHandle> side_verts;
1044 11904 : for (const auto i : make_range(side_elem->n_nodes()))
1045 : {
1046 8928 : const auto node_it = _node_id_to_handle.find(side_elem->node_id(i));
1047 : mooseAssert(node_it != _node_id_to_handle.end(),
1048 : "Could not map libMesh node " + std::to_string(side_elem->node_id(i)) +
1049 : " to a MOAB vertex while assigning DAGMC boundary conditions.");
1050 8928 : side_verts.insert(node_it->second);
1051 : }
1052 :
1053 : // Retrieve the MOAB triangle faces adjacent to this tet
1054 : moab::Range tri_faces;
1055 2976 : check(_moab->get_adjacencies(&tet, num_entities, surface_dimension, false, tri_faces));
1056 :
1057 10676 : for (const auto tri : tri_faces)
1058 : {
1059 3850 : if (skin.find(tri) == skin.end())
1060 1930 : continue;
1061 :
1062 2588 : const moab::EntityHandle * conn = nullptr;
1063 2588 : int nconn = 0;
1064 2588 : check(_moab->get_connectivity(tri, conn, nconn));
1065 :
1066 : // A skin triangle lies on this libMesh side if all three of its MOAB vertices
1067 : // are contained in the side's vertex set.
1068 2588 : if (!side_verts.count(conn[0]) || !side_verts.count(conn[1]) || !side_verts.count(conn[2]))
1069 668 : continue;
1070 :
1071 1920 : switch (bc_type)
1072 : {
1073 : case BoundaryConditionType::Transmission:
1074 : // Nothing to do - tris stay in transmission_tris
1075 : break;
1076 448 : case BoundaryConditionType::Vacuum:
1077 448 : transmission_tris.erase(tri);
1078 448 : vacuum_tris.insert(tri);
1079 : break;
1080 1472 : case BoundaryConditionType::Reflective:
1081 1472 : transmission_tris.erase(tri);
1082 1472 : reflective_tris.insert(tri);
1083 : break;
1084 0 : default:
1085 0 : mooseError("Unhandled boundary condition type!");
1086 : }
1087 : }
1088 2976 : }
1089 : }
1090 : }
1091 :
1092 : void
1093 3693 : MoabSkinner::createSurfacesFromSkin(const moab::Range & region,
1094 : moab::Range & skin,
1095 : VolData & voldata,
1096 : unsigned int & surf_id)
1097 : {
1098 : moab::Range transmission_tris, vacuum_tris, reflective_tris;
1099 3693 : splitSkinByBoundaryCondition(region, skin, transmission_tris, vacuum_tris, reflective_tris);
1100 :
1101 : // Create surfaces for each BC class separately so every surface meshset carries
1102 : // one BC type. Transmission is the DAGMC default and needs no group.
1103 3692 : createSurfaces(transmission_tris, voldata, surf_id, BoundaryConditionType::Transmission);
1104 3692 : createSurfaces(vacuum_tris, voldata, surf_id, BoundaryConditionType::Vacuum);
1105 3692 : createSurfaces(reflective_tris, voldata, surf_id, BoundaryConditionType::Reflective);
1106 3691 : }
1107 :
1108 : void
1109 5235 : MoabSkinner::recordBoundaryConditionSurface(moab::EntityHandle surface_set,
1110 : BoundaryConditionType bc_type)
1111 : {
1112 : // transmission is the DAGMC default - we don't need to tag the transmission surface
1113 5235 : if (bc_type == BoundaryConditionType::Transmission)
1114 5207 : return;
1115 :
1116 28 : const auto [it, inserted] = _surface_bc_types.emplace(surface_set, bc_type);
1117 28 : if (!inserted && it->second != bc_type)
1118 1 : mooseError("A DAGMC surface was assigned two different boundary conditions ('",
1119 1 : boundaryConditionGroupName(it->second),
1120 : "' and '",
1121 1 : boundaryConditionGroupName(bc_type),
1122 : "').");
1123 : }
1124 :
1125 : MoabSkinner::BoundaryConditionType
1126 2414 : MoabSkinner::recordedBoundaryCondition(moab::EntityHandle surface_set) const
1127 : {
1128 : const auto it = _surface_bc_types.find(surface_set);
1129 2414 : return it == _surface_bc_types.end() ? BoundaryConditionType::Transmission : it->second;
1130 : }
1131 :
1132 : std::string
1133 16 : MoabSkinner::boundaryConditionGroupName(BoundaryConditionType bc_type) const
1134 : {
1135 16 : switch (bc_type)
1136 : {
1137 : case BoundaryConditionType::Vacuum:
1138 8 : return "boundary:Vacuum";
1139 : case BoundaryConditionType::Reflective:
1140 8 : return "boundary:Reflecting";
1141 : case BoundaryConditionType::Transmission:
1142 : mooseAssert(bc_type != BoundaryConditionType::Transmission,
1143 : "Surface set groups should not be created for transmission BCs since "
1144 : "transmission is the DAGMC default.");
1145 0 : return "";
1146 0 : default:
1147 0 : mooseError("No DAGMC boundary condition group exists for this boundary condition type");
1148 : }
1149 : }
1150 :
1151 : unsigned int
1152 8 : MoabSkinner::firstBoundaryConditionGroupID() const
1153 : {
1154 : // IDs used by material groups.
1155 : // buildGraveyard() and the implicit complement each consume one additional ID when enabled.
1156 8 : unsigned int gid = nBins() + 1;
1157 8 : gid += _build_graveyard + _set_implicit_complement_material;
1158 8 : return gid;
1159 : }
1160 :
1161 : void
1162 70 : MoabSkinner::createBoundaryConditionGroups()
1163 : {
1164 70 : if (!_set_bcs)
1165 60 : return;
1166 :
1167 : std::map<BoundaryConditionType, std::vector<moab::EntityHandle>> surfaces_by_type;
1168 36 : for (const auto & [surf, bc_type] : _surface_bc_types)
1169 26 : surfaces_by_type[bc_type].push_back(surf);
1170 :
1171 17 : if (!_vacuum_bcs_surface_ids.empty() && !surfaces_by_type.count(BoundaryConditionType::Vacuum))
1172 1 : paramError("vacuum_bcs_surfaces",
1173 : "'vacuum_bcs_surfaces' was specified but no skinned DAGMC surfaces were assigned "
1174 : "vacuum boundary conditions. Verify the sideset names or IDs correspond to "
1175 : "boundary faces of the mesh.");
1176 :
1177 9 : if (!_reflective_bcs_surface_ids.empty() &&
1178 8 : !surfaces_by_type.count(BoundaryConditionType::Reflective))
1179 1 : paramError("reflective_bcs_surfaces",
1180 : "'reflective_bcs_surfaces' was specified but no skinned DAGMC surfaces were "
1181 : "assigned reflective boundary conditions. Verify the sideset names or IDs "
1182 : "correspond to boundary faces of the mesh.");
1183 :
1184 8 : unsigned int gid = firstBoundaryConditionGroupID();
1185 22 : for (const auto & [bc_type, surfs] : surfaces_by_type)
1186 : {
1187 14 : moab::EntityHandle group = 0;
1188 14 : createGroup(gid++, boundaryConditionGroupName(bc_type), group);
1189 40 : for (const auto surf_set : surfs)
1190 26 : check(_moab->add_entities(group, &surf_set, 1));
1191 : }
1192 : }
1193 :
1194 : void
1195 72 : MoabSkinner::findSurfaces()
1196 : {
1197 : _surface_bc_types.clear();
1198 :
1199 : // Find all neighbours in mesh
1200 72 : getDAGMCGeometryMesh().find_neighbors();
1201 :
1202 : // Counter for volumes
1203 72 : unsigned int vol_id = 0;
1204 :
1205 : // Counter for surfaces
1206 72 : unsigned int surf_id = 0;
1207 :
1208 : // Loop over material bins
1209 301 : for (unsigned int iMat = 0; iMat < _n_block_bins; iMat++)
1210 : {
1211 : // Loop over density bins
1212 630 : for (unsigned int iDen = 0; iDen < _n_density_bins; iDen++)
1213 : {
1214 : // Loop over temperature bins
1215 1456 : for (unsigned int iVar = 0; iVar < _n_temperature_bins; iVar++)
1216 : {
1217 : // Update material name
1218 1057 : auto updated_mat_name = materialName(iMat, iDen, iVar);
1219 :
1220 : // Create a material group
1221 1057 : int iSortBin = getBin(iVar, iDen, iMat);
1222 :
1223 : // For DagMC to fill a cell with a material, we first create a group
1224 : // with that name, and then assign it with createVol (called inside findSurface)
1225 : moab::EntityHandle group_set;
1226 1057 : unsigned int group_id = iSortBin + 1;
1227 1057 : createGroup(group_id, updated_mat_name, group_set);
1228 :
1229 : // Sort elems in this mat-density-temp bin into local regions
1230 : std::vector<moab::Range> regions;
1231 3171 : groupLocalElems(_elem_bins.at(iSortBin), regions);
1232 :
1233 : // Loop over all regions and find surfaces
1234 2902 : for (const auto & region : regions)
1235 : {
1236 : moab::EntityHandle volume_set;
1237 1847 : findSurface(region, group_set, vol_id, surf_id, volume_set);
1238 : }
1239 1055 : }
1240 : }
1241 : }
1242 :
1243 70 : if (_build_graveyard)
1244 48 : buildGraveyard(vol_id, surf_id);
1245 :
1246 70 : if (_set_implicit_complement_material)
1247 : {
1248 : moab::EntityHandle comp_group;
1249 2 : unsigned int comp_id = nBins() + 1 + _build_graveyard;
1250 2 : createGroup(comp_id, _implicit_complement_group_name, comp_group);
1251 2 : moab::EntityHandle arbitray_volume = 0;
1252 2 : for (const auto & surf_pair : surfsToVols)
1253 : {
1254 : const auto & vols = surf_pair.second;
1255 2 : arbitray_volume = vols.front().vol;
1256 2 : break;
1257 : }
1258 2 : check(_moab->add_entities(comp_group, &arbitray_volume, 1));
1259 : }
1260 :
1261 70 : createBoundaryConditionGroups();
1262 :
1263 : // Write MOAB volume and/or skin meshes to file
1264 68 : write();
1265 68 : }
1266 :
1267 : void
1268 68 : MoabSkinner::write()
1269 : {
1270 : // Only write to file on root process
1271 68 : if (processor_id() != 0)
1272 0 : return;
1273 :
1274 68 : std::string extension = std::to_string(_n_write) + ".h5m";
1275 :
1276 68 : if (_output_skins)
1277 : {
1278 : // Generate list of surfaces to write
1279 : std::vector<moab::EntityHandle> surfs;
1280 3549 : for (const auto & itsurf : surfsToVols)
1281 3532 : surfs.push_back(itsurf.first);
1282 :
1283 17 : std::string filename = "moab_skins_" + extension;
1284 :
1285 17 : if (_verbose)
1286 17 : _console << "Writing MOAB skins to " << filename << "...";
1287 :
1288 17 : check(_moab->write_mesh(filename.c_str(), surfs.data(), surfs.size()));
1289 17 : }
1290 :
1291 68 : if (_output_full)
1292 : {
1293 2 : std::string filename = "moab_mesh_" + extension;
1294 :
1295 2 : if (_verbose)
1296 1 : _console << "Writing MOAB mesh to " << filename << std::endl;
1297 :
1298 2 : check(_moab->write_mesh(filename.c_str()));
1299 : }
1300 :
1301 68 : _n_write++;
1302 : }
1303 :
1304 : void
1305 1057 : MoabSkinner::groupLocalElems(std::set<dof_id_type> elems, std::vector<moab::Range> & localElems)
1306 : {
1307 2904 : while (!elems.empty())
1308 : {
1309 :
1310 : // Create a new local range of moab handles
1311 : moab::Range local;
1312 :
1313 : // Retrieve and remove the fisrt elem
1314 : auto it = elems.begin();
1315 1847 : dof_id_type next = *it;
1316 1847 : elems.erase(it);
1317 :
1318 : std::set<dof_id_type> neighbors;
1319 1847 : neighbors.insert(next);
1320 :
1321 11657 : while (!neighbors.empty())
1322 : {
1323 :
1324 : std::set<dof_id_type> new_neighbors;
1325 :
1326 : // Loop over all the new neighbors
1327 93443 : for (auto & next : neighbors)
1328 : {
1329 :
1330 : // Get the MOAB handles, and add to local set
1331 : // (May be more than one if this libMesh elem has sub-tetrahedra)
1332 83633 : if (_id_to_elem_handles.find(next) == _id_to_elem_handles.end())
1333 0 : mooseError("No entity handles found for libmesh id.");
1334 :
1335 83633 : std::vector<moab::EntityHandle> ents = _id_to_elem_handles[next];
1336 178578 : for (const auto ent : ents)
1337 94945 : local.insert(ent);
1338 :
1339 : // Get the libMesh element
1340 83633 : Elem & elem = getDAGMCGeometryMesh().elem_ref(next);
1341 :
1342 : // How many nearest neighbors (general element)?
1343 : unsigned int NN = elem.n_neighbors();
1344 :
1345 : // Loop over neighbors
1346 418165 : for (unsigned int i = 0; i < NN; i++)
1347 : {
1348 : const Elem * nnptr = elem.neighbor_ptr(i);
1349 : // If on boundary, some may be null ptrs
1350 334532 : if (nnptr == nullptr)
1351 24500 : continue;
1352 :
1353 310032 : dof_id_type idnn = nnptr->id();
1354 :
1355 : // Select only those that are in the current bin
1356 310032 : if (elems.find(idnn) != elems.end())
1357 : {
1358 81786 : new_neighbors.insert(idnn);
1359 : // Remove from those still available
1360 : elems.erase(idnn);
1361 : }
1362 : }
1363 83633 : }
1364 :
1365 : // Found all the new neighbors, done with current set.
1366 : neighbors = new_neighbors;
1367 : }
1368 :
1369 : // Save this moab range of local neighbors
1370 1847 : localElems.push_back(local);
1371 : }
1372 1057 : }
1373 :
1374 : void
1375 76 : MoabSkinner::reset()
1376 : {
1377 76 : _moab.reset(new moab::Core());
1378 76 : skinner.reset(new moab::Skinner(_moab.get()));
1379 76 : gtt.reset(new moab::GeomTopoTool(_moab.get()));
1380 :
1381 : _tet_mesh.reset();
1382 76 : _tet_mesh_built = false;
1383 :
1384 : // Clear entity set maps
1385 : surfsToVols.clear();
1386 76 : }
1387 :
1388 : unsigned int
1389 80903 : MoabSkinner::getBin(const unsigned int & i_temp,
1390 : const unsigned int & i_density,
1391 : const unsigned int & i_block) const
1392 : {
1393 80903 : return _n_temperature_bins * (_n_density_bins * i_block + i_density) + i_temp;
1394 : }
1395 :
1396 : void
1397 1847 : MoabSkinner::findSurface(const moab::Range & region,
1398 : moab::EntityHandle group,
1399 : unsigned int & vol_id,
1400 : unsigned int & surf_id,
1401 : moab::EntityHandle & volume_set)
1402 : {
1403 : // Create a volume set
1404 1847 : vol_id++;
1405 1847 : createVol(vol_id, volume_set, group);
1406 :
1407 : // Find surfaces from these regions
1408 : moab::Range tris; // The tris of the surfaces
1409 : moab::Range rtris; // The tris which are reversed with respect to their surfaces
1410 1847 : skinner->find_skin(0, region, false, tris, &rtris);
1411 :
1412 : // Create surface sets, classifying by boundary condition. BC sorting happens here,
1413 : // while the current region and its skin result are in hand, rather than in a
1414 : // separate post-processing pass.
1415 1847 : VolData vdata = {volume_set, Sense::FORWARDS};
1416 1847 : createSurfacesFromSkin(region, tris, vdata, surf_id);
1417 :
1418 : // Create surface sets for the reversed tris
1419 1846 : vdata.sense = Sense::BACKWARDS;
1420 1846 : createSurfacesFromSkin(region, rtris, vdata, surf_id);
1421 1845 : }
1422 :
1423 : void
1424 11076 : MoabSkinner::createSurfaces(moab::Range & faces,
1425 : VolData & voldata,
1426 : unsigned int & surf_id,
1427 : BoundaryConditionType bc_type)
1428 : {
1429 11076 : if (faces.empty())
1430 : return;
1431 :
1432 : // Loop over the surfaces we have already created
1433 1800979 : for (const auto & surfpair : surfsToVols)
1434 : {
1435 : // Local copies of surf/vols
1436 1799118 : moab::EntityHandle surf = surfpair.first;
1437 1799118 : std::vector<VolData> vols = surfpair.second;
1438 :
1439 : // First get the entities in this surface
1440 : moab::Range tris;
1441 1799118 : check(_moab->get_entities_by_handle(surf, tris));
1442 :
1443 : // Find any tris that live in both surfs
1444 1799118 : moab::Range overlap = moab::intersect(tris, faces);
1445 1799118 : if (!overlap.empty())
1446 : {
1447 : // Check if the tris are a subset or the entire surf
1448 3278 : if (tris.size() == overlap.size())
1449 : {
1450 : // Whole surface -> just update the volume relationships and BC record
1451 864 : updateSurfData(surf, voldata);
1452 864 : recordBoundaryConditionSurface(surf, bc_type);
1453 : }
1454 : else
1455 : {
1456 : // Overlap is a subset: remove shared tris from this surface and create a new
1457 : // shared surface carrying both volume relationships and the BC type
1458 2414 : check(_moab->remove_entities(surf, overlap));
1459 :
1460 : // Append our new volume to the list that share this surf
1461 2414 : vols.push_back(voldata);
1462 :
1463 : // The shared tris may have been assigned a BC when 'surf' was created (e.g. by
1464 : // the region on the other side of an internal surface); merge that record with
1465 : // the current classification so the BC is not lost when the tris move to the
1466 : // new shared surface
1467 : auto merged_bc = bc_type;
1468 2414 : const auto existing_bc = recordedBoundaryCondition(surf);
1469 2414 : if (merged_bc == BoundaryConditionType::Transmission)
1470 : merged_bc = existing_bc;
1471 0 : else if (existing_bc != BoundaryConditionType::Transmission && existing_bc != merged_bc)
1472 0 : mooseError("A DAGMC surface was assigned both vacuum and reflective boundary "
1473 : "conditions. This surface is shared between two skinned regions "
1474 : "(e.g. an internal surface between two blocks) and received a "
1475 : "different boundary condition from each side. Check 'vacuum_bcs_surfaces' "
1476 : "and 'reflective_bcs_surfaces' for sidesets that cover the same mesh faces "
1477 : "from opposite sides.");
1478 :
1479 : moab::EntityHandle shared_surf;
1480 2414 : surf_id++;
1481 2414 : createSurf(surf_id, shared_surf, overlap, vols, merged_bc);
1482 : }
1483 :
1484 : // Subtract from the input list
1485 24149 : for (auto & shared : overlap)
1486 20872 : faces.erase(shared);
1487 :
1488 3277 : if (faces.empty())
1489 : break;
1490 : }
1491 1799117 : }
1492 :
1493 3590 : if (!faces.empty())
1494 : {
1495 : moab::EntityHandle surface_set;
1496 1861 : std::vector<VolData> voldatavec(1, voldata);
1497 1861 : surf_id++;
1498 1861 : createSurf(surf_id, surface_set, faces, voldatavec, bc_type);
1499 1861 : }
1500 : }
1501 :
1502 : void
1503 48 : MoabSkinner::buildGraveyard(unsigned int & vol_id, unsigned int & surf_id)
1504 : {
1505 : // Create the graveyard group
1506 : moab::EntityHandle graveyard;
1507 48 : unsigned int id = nBins() + 1;
1508 48 : createGroup(id, "mat:Graveyard", graveyard);
1509 :
1510 : // Create a volume set
1511 : moab::EntityHandle volume_set;
1512 48 : createVol(++vol_id, volume_set, graveyard);
1513 :
1514 : // Set up for the volume data to pass to surfs
1515 48 : VolData vdata = {volume_set, Sense::FORWARDS};
1516 :
1517 : // Find a bounding box
1518 48 : BoundingBox bbox = MeshTools::create_bounding_box(getMooseMesh().getMesh());
1519 :
1520 : // Build the two cubic surfaces defining the graveyard
1521 48 : createSurfaceFromBox(
1522 : bbox, vdata, surf_id, false /* normals point into box */, _graveyard_scale_inner);
1523 48 : createSurfaceFromBox(
1524 : bbox, vdata, surf_id, true /* normals point out of box */, _graveyard_scale_outer);
1525 48 : }
1526 :
1527 : void
1528 96 : MoabSkinner::createSurfaceFromBox(const BoundingBox & box,
1529 : const VolData & voldata,
1530 : unsigned int & surf_id,
1531 : bool normalout,
1532 : const Real & factor)
1533 : {
1534 96 : std::vector<moab::EntityHandle> vert_handles = createNodesFromBox(box, factor);
1535 :
1536 : // Create the tris in 4 groups of 3 (4 open tetrahedra)
1537 : moab::Range tris;
1538 96 : createCornerTris(vert_handles, 0, 1, 2, 4, normalout, tris);
1539 96 : createCornerTris(vert_handles, 3, 2, 1, 7, normalout, tris);
1540 96 : createCornerTris(vert_handles, 6, 4, 2, 7, normalout, tris);
1541 96 : createCornerTris(vert_handles, 5, 1, 4, 7, normalout, tris);
1542 :
1543 : moab::EntityHandle surface_set;
1544 96 : std::vector<VolData> voldatavec(1, voldata);
1545 96 : surf_id++;
1546 96 : createSurf(surf_id, surface_set, tris, voldatavec);
1547 192 : }
1548 :
1549 : std::vector<moab::EntityHandle>
1550 96 : MoabSkinner::createNodesFromBox(const BoundingBox & box, const Real & factor) const
1551 : {
1552 : std::vector<moab::EntityHandle> vert_handles;
1553 :
1554 : // Fetch the vertices of the box
1555 96 : auto verts = geom_utils::boxCorners(box, factor);
1556 :
1557 : // Array to represent a coord in MOAB
1558 : double coord[3];
1559 :
1560 : // Create the vertices in MOAB and get the handles
1561 864 : for (const auto & vert : verts)
1562 : {
1563 768 : coord[0] = vert(0) * _scaling;
1564 768 : coord[1] = vert(1) * _scaling;
1565 768 : coord[2] = vert(2) * _scaling;
1566 :
1567 : moab::EntityHandle ent;
1568 768 : check(_moab->create_vertex(coord, ent));
1569 768 : vert_handles.push_back(ent);
1570 : }
1571 :
1572 96 : return vert_handles;
1573 96 : }
1574 :
1575 : void
1576 384 : MoabSkinner::createCornerTris(const std::vector<moab::EntityHandle> & verts,
1577 : unsigned int corner,
1578 : unsigned int v1,
1579 : unsigned int v2,
1580 : unsigned int v3,
1581 : bool normalout,
1582 : moab::Range & surface_tris)
1583 : {
1584 : // Create 3 tris stemming from one corner (i.e. an open tetrahedron)
1585 : // Assume first is the central corner, and the others are labelled clockwise looking down on the
1586 : // corner
1587 384 : unsigned int indices[3] = {v1, v2, v3};
1588 :
1589 : // Create each tri by a cyclic permutation of indices
1590 : // Values of i1, i2 in the loop: 0,1; 1,2; 2;0
1591 1536 : for (unsigned int i = 0; i < 3; i++)
1592 : {
1593 1152 : int i1 = indices[i % 3];
1594 1152 : int i2 = indices[(i + 1) % 3];
1595 1152 : if (normalout) // anti-clockwise: normal points outwards
1596 576 : surface_tris.insert(createTri(verts, corner, i2, i1));
1597 : else // clockwise: normal points inwards
1598 576 : surface_tris.insert(createTri(verts, corner, i1, i2));
1599 : }
1600 384 : }
1601 :
1602 : moab::EntityHandle
1603 1152 : MoabSkinner::createTri(const std::vector<moab::EntityHandle> & vertices,
1604 : unsigned int v1,
1605 : unsigned int v2,
1606 : unsigned int v3)
1607 : {
1608 : moab::EntityHandle triangle;
1609 1152 : moab::EntityHandle connectivity[3] = {vertices[v1], vertices[v2], vertices[v3]};
1610 1152 : check(_moab->create_element(moab::MBTRI, connectivity, 3, triangle));
1611 1152 : return triangle;
1612 : }
1613 :
1614 : void
1615 37 : MoabSkinner::setGraveyard(bool build)
1616 : {
1617 37 : if (build != _build_graveyard)
1618 : {
1619 34 : std::string original = _build_graveyard ? "true" : "false";
1620 34 : std::string change = _build_graveyard ? "false" : "true";
1621 17 : mooseWarning("Overriding graveyard setting from ",
1622 : original,
1623 : " to ",
1624 : change,
1625 : ".\n"
1626 : "To hide this warning, set 'build_graveyard = ",
1627 : change,
1628 : "'");
1629 : }
1630 :
1631 36 : _build_graveyard = build;
1632 36 : }
1633 : #endif
|
