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 :
11 : #include "libmesh/elem.h"
12 : #include "libmesh/enum_io_package.h"
13 : #include "libmesh/enum_order.h"
14 : #include "libmesh/enum_fe_family.h"
15 : #include "libmesh/equation_systems.h"
16 : #include "libmesh/system.h"
17 : #include "libmesh/mesh_tools.h"
18 :
19 : registerMooseObject("CardinalApp", MoabSkinner);
20 :
21 : InputParameters
22 88 : MoabSkinner::validParams()
23 : {
24 88 : InputParameters params = GeneralUserObject::validParams();
25 176 : params.addParam<bool>("verbose", false, "Whether to print diagnostic information");
26 :
27 : // temperature binning
28 176 : params.addRequiredParam<std::string>("temperature",
29 : "Temperature variable by which to bin elements");
30 264 : params.addRangeCheckedParam<Real>(
31 176 : "temperature_min", 0.0, "temperature_min >= 0.0", "Lower bound of temperature bins");
32 176 : params.addRequiredParam<Real>("temperature_max", "Upper bound of temperature bins");
33 176 : params.addRequiredRangeCheckedParam<unsigned int>(
34 : "n_temperature_bins", "n_temperature_bins > 0", "Number of temperature bins");
35 :
36 : // density binning
37 176 : params.addParam<std::string>("density", "Density variable by which to bin elements");
38 264 : params.addRangeCheckedParam<Real>(
39 176 : "density_min", 0.0, "density_min >= 0.0", "Lower bound of density bins");
40 176 : params.addParam<Real>("density_max", "Upper bound of density bins");
41 176 : params.addRangeCheckedParam<unsigned int>(
42 : "n_density_bins", "n_density_bins > 0", "Number of density bins");
43 :
44 176 : params.addParam<std::vector<std::string>>(
45 : "material_names",
46 : "List of names for each subdomain to use for naming the new volumes created in MOAB. "
47 : "You only need to set this if using this skinner independent of OpenMC; otherwise, "
48 : "these names are auto-deduced from OpenMC");
49 :
50 264 : params.addRangeCheckedParam<Real>(
51 176 : "faceting_tol", 1e-4, "faceting_tol > 0", "Faceting tolerance for DagMC");
52 264 : params.addRangeCheckedParam<Real>(
53 176 : "geom_tol", 1e-6, "geom_tol > 0", "Geometry tolerance for DagMC");
54 :
55 176 : params.addParam<bool>(
56 176 : "build_graveyard", false, "Whether to build a graveyard around the geometry");
57 264 : params.addRangeCheckedParam<Real>(
58 : "graveyard_scale_inner",
59 176 : 1.01,
60 : "graveyard_scale_inner > 1",
61 : "Multiplier on mesh bounding box to form inner graveyard surface");
62 176 : params.addParam<Real>("graveyard_scale_outer",
63 176 : 1.10,
64 : "Multiplier on mesh bounding box to form outer graveyard surface");
65 176 : params.addParam<std::string>("implicit_complement_material",
66 : "Assigns OpenMC material name or ID to the implicit complement "
67 : "region. If not provided, void material is assigned by default.");
68 : // TODO: would be nice to support other file formats as well, like exodus
69 176 : params.addParam<bool>(
70 : "output_skins",
71 176 : false,
72 : "Whether the skinned MOAB mesh (skins generated from the "
73 : "libMesh [Mesh]) should be written to a file. The files will be named moab_skins_<n>.h5m, "
74 : "where <n> is the time step index. You can then visualize these files by running "
75 : "'mbconvert'.");
76 176 : params.addParam<bool>("output_full",
77 176 : false,
78 : "Whether the MOAB mesh (copied from the libMesh [Mesh]) should "
79 : "be written to a file. The files will be named moab_full_<n>.h5m, where "
80 : "<n> is the time step index. "
81 : "You can then visualize these files by running 'mbconvert'.");
82 176 : params.addParam<bool>("use_displaced_mesh",
83 176 : false,
84 : "Whether the skinned mesh should be generated from a displaced mesh ");
85 88 : params.addClassDescription("Re-generate the OpenMC geometry on-the-fly according to changes in "
86 : "the mesh geometry and/or contours in temperature and density");
87 88 : return params;
88 0 : }
89 :
90 44 : MoabSkinner::MoabSkinner(const InputParameters & parameters)
91 : : GeneralUserObject(parameters),
92 44 : _serialized_solution(NumericVector<Number>::build(_communicator).release()),
93 88 : _verbose(getParam<bool>("verbose")),
94 88 : _temperature_name(getParam<std::string>("temperature")),
95 88 : _temperature_min(getParam<Real>("temperature_min")),
96 88 : _temperature_max(getParam<Real>("temperature_max")),
97 88 : _n_temperature_bins(getParam<unsigned int>("n_temperature_bins")),
98 44 : _temperature_bin_width((_temperature_max - _temperature_min) / _n_temperature_bins),
99 88 : _bin_by_density(isParamValid("density")),
100 44 : _faceting_tol(getParam<Real>("faceting_tol")),
101 88 : _geom_tol(getParam<Real>("geom_tol")),
102 88 : _graveyard_scale_inner(getParam<double>("graveyard_scale_inner")),
103 88 : _graveyard_scale_outer(getParam<double>("graveyard_scale_outer")),
104 88 : _output_skins(getParam<bool>("output_skins")),
105 88 : _output_full(getParam<bool>("output_full")),
106 44 : _scaling(1.0),
107 44 : _n_write(0),
108 88 : _standalone(true)
109 : {
110 88 : _build_graveyard = getParam<bool>("build_graveyard");
111 88 : _use_displaced = getParam<bool>("use_displaced_mesh");
112 :
113 88 : if (isParamSetByUser("implicit_complement_material"))
114 : {
115 : // If the user specify a material that doesn't exist in materials.xml file, OpenMC
116 : // will catch the mistake.
117 2 : _set_implicit_complement_material = true;
118 : _implicit_complement_group_name =
119 8 : "mat:" + getParam<std::string>("implicit_complement_material") + "_comp";
120 : }
121 :
122 : // we can probably support this in the future, it's just not implemented yet
123 44 : if (!getMooseMesh().getMesh().is_serial())
124 0 : mooseError("MoabSkinner does not yet support distributed meshes!");
125 :
126 : // Create MOAB interface
127 44 : _moab = std::make_shared<moab::Core>();
128 :
129 : // Create a skinner
130 88 : skinner = std::make_unique<moab::Skinner>(_moab.get());
131 :
132 : // Create a geom topo tool
133 88 : gtt = std::make_unique<moab::GeomTopoTool>(_moab.get());
134 :
135 44 : if (_bin_by_density)
136 : {
137 40 : checkRequiredParam(parameters, "density_min", "binning by density");
138 40 : checkRequiredParam(parameters, "density_max", "binning by density");
139 40 : checkRequiredParam(parameters, "n_density_bins", "binning by density");
140 :
141 40 : _density_min = getParam<Real>("density_min");
142 40 : _density_max = getParam<Real>("density_max");
143 40 : _n_density_bins = getParam<unsigned int>("n_density_bins");
144 40 : _density_name = getParam<std::string>("density");
145 20 : _density_bin_width = (_density_max - _density_min) / _n_density_bins;
146 :
147 20 : if (_density_max < _density_min)
148 1 : paramError("density_max", "'density_max' must be greater than 'density_min'");
149 : }
150 : else
151 : {
152 48 : checkUnusedParam(parameters, "density_min", "not binning by density");
153 48 : checkUnusedParam(parameters, "density_max", "not binning by density");
154 48 : checkUnusedParam(parameters, "n_density_bins", "not binning by density");
155 :
156 24 : _n_density_bins = 1;
157 : }
158 :
159 43 : if (_build_graveyard)
160 : {
161 32 : if (_graveyard_scale_outer < _graveyard_scale_inner)
162 1 : paramError("graveyard_scale_outer",
163 : "'graveyard_scale_outer' must be greater than 'graveyard_scale_inner'!");
164 : }
165 : else
166 : {
167 22 : checkUnusedParam(parameters, "graveyard_scale_inner", "'build_graveyard' is false");
168 22 : checkUnusedParam(parameters, "graveyard_scale_outer", "'build_graveyard' is false");
169 : }
170 :
171 : // get variable numbers
172 42 : _temperature_var_num = getAuxiliaryVariableNumber(_temperature_name, "temperature");
173 40 : if (_bin_by_density)
174 : {
175 17 : if (_temperature_name == _density_name)
176 1 : mooseError("The 'temperature' and 'density' variables cannot be the same!");
177 :
178 16 : _density_var_num = getAuxiliaryVariableNumber(_density_name, "density");
179 : }
180 :
181 38 : if (_temperature_max <= _temperature_min)
182 1 : paramError("temperature_max", "'temperature_max' must be greater than 'temperature_min'");
183 :
184 205 : for (unsigned int i = 0; i < _n_temperature_bins + 1; ++i)
185 168 : _temperature_bin_bounds.push_back(_temperature_min + i * _temperature_bin_width);
186 :
187 153 : for (unsigned int i = 0; i < _n_density_bins + 1; ++i)
188 116 : _density_bin_bounds.push_back(_density_min + i * _density_bin_width);
189 :
190 : // node numberings for first-order tets
191 74 : _tet4_nodes.push_back({0, 1, 2, 3});
192 :
193 : // node numbers for second-order tets
194 74 : _tet10_nodes.push_back({0, 4, 6, 7});
195 74 : _tet10_nodes.push_back({1, 5, 4, 8});
196 74 : _tet10_nodes.push_back({2, 6, 5, 9});
197 74 : _tet10_nodes.push_back({7, 8, 9, 3});
198 74 : _tet10_nodes.push_back({4, 9, 7, 8});
199 74 : _tet10_nodes.push_back({4, 5, 9, 8});
200 74 : _tet10_nodes.push_back({4, 7, 9, 6});
201 74 : _tet10_nodes.push_back({4, 9, 5, 6});
202 :
203 37 : moab::MBErrorHandler_Init();
204 37 : }
205 :
206 : void
207 42 : MoabSkinner::finalize()
208 : {
209 42 : moab::MBErrorHandler_Finalize();
210 42 : }
211 :
212 : moab::ErrorCode
213 1949932 : MoabSkinner::check(const moab::ErrorCode input) const
214 : {
215 : #ifdef DEBUG
216 : MB_CHK_ERR(input);
217 : #endif
218 1949932 : return moab::MB_SUCCESS;
219 : }
220 :
221 : unsigned int
222 58 : MoabSkinner::getAuxiliaryVariableNumber(const std::string & name,
223 : const std::string & param_name) const
224 : {
225 58 : if (!_fe_problem.getAuxiliarySystem().hasVariable(name))
226 2 : paramError(param_name, "Cannot find auxiliary variable '", name, "'!");
227 :
228 : // we require these variables to be constant monomial
229 56 : auto type = _fe_problem.getAuxiliarySystem().getFieldVariable<Real>(0, name).feType();
230 56 : if (type.family != MONOMIAL || type.order != 0)
231 1 : paramError(param_name, "Auxiliary variable '", name, "' must be a CONSTANT MONOMIAL type!");
232 :
233 55 : return _fe_problem.getAuxiliarySystem().getFieldVariable<Real>(0, name).number();
234 : }
235 :
236 : MooseMesh &
237 77987 : MoabSkinner::getMooseMesh()
238 : {
239 103790 : if (_use_displaced && _fe_problem.getDisplacedProblem() == nullptr)
240 0 : mooseError("Displaced mesh was requested but the displaced problem does not exist. "
241 : "set use_displaced_mesh = False");
242 25803 : return ((_use_displaced && _fe_problem.getDisplacedProblem())
243 181777 : ? _fe_problem.getDisplacedProblem()->mesh()
244 130171 : : _fe_problem.mesh());
245 : }
246 :
247 : void
248 112 : MoabSkinner::initialize()
249 : {
250 112 : findBlocks();
251 :
252 112 : if (_standalone)
253 : {
254 60 : checkRequiredParam(
255 : parameters(), "material_names", "using skinner independent of an OpenMC [Problem]");
256 90 : _material_names = getParam<std::vector<std::string>>("material_names");
257 :
258 30 : if (_material_names.size() != _n_block_bins)
259 2 : paramError("material_names",
260 : "This parameter must be the same length as the number of "
261 1 : "subdomains in the mesh (" +
262 : Moose::stringify(_n_block_bins) + ")");
263 : }
264 : else
265 : {
266 : // the external class should have set the value of _material_names
267 164 : checkUnusedParam(parameters(),
268 : "material_names",
269 : "using the skinner in conjunction with an OpenMC [Problem]");
270 : }
271 :
272 : // Set spatial dimension in MOAB
273 111 : check(_moab->set_dimension(getMooseMesh().getMesh().spatial_dimension()));
274 :
275 : // Create a meshset representing all of the MOAB tets
276 111 : check(_moab->create_meshset(moab::MESHSET_SET, _all_tets));
277 :
278 111 : createTags();
279 :
280 111 : createMOABElems();
281 110 : }
282 :
283 : void
284 46 : MoabSkinner::execute()
285 : {
286 46 : if (_standalone)
287 14 : update();
288 42 : }
289 :
290 : void
291 19 : MoabSkinner::setUseDisplacedMesh(const bool & use)
292 : {
293 19 : if ((use != _use_displaced) && isParamSetByUser("use_displaced_mesh"))
294 0 : mooseWarning("Overriding 'use_displaced_mesh' to " + std::to_string(use) +
295 : " to match the displaced problem action.");
296 19 : _use_displaced = use;
297 19 : }
298 :
299 : void
300 46 : MoabSkinner::update()
301 : {
302 46 : _console << "Skinning geometry into " << _n_temperature_bins << " temperature bins, "
303 46 : << _n_density_bins << " density bins, and " << _n_block_bins << " block bins... "
304 46 : << std::endl;
305 :
306 46 : if (_use_displaced && _standalone)
307 : {
308 : // we are responsible for updating the mesh if running in standalone mode; otherwise, the
309 : // OpenMCCellAverageProblem class does it
310 4 : _fe_problem.getDisplacedProblem()->updateMesh();
311 : }
312 :
313 : // Clear MOAB mesh data from last timestep
314 46 : reset();
315 :
316 46 : _serialized_solution->init(_fe_problem.getAuxiliarySystem().sys().n_dofs(), false, SERIAL);
317 46 : _fe_problem.getAuxiliarySystem().solution().localize(*_serialized_solution);
318 :
319 : // Re-initialise the mesh data
320 46 : initialize();
321 :
322 : // Sort libMesh elements into bins
323 46 : sortElemsByResults();
324 :
325 : // Find the surfaces of local temperature regions
326 42 : findSurfaces();
327 42 : }
328 :
329 : void
330 112 : MoabSkinner::findBlocks()
331 : {
332 : _blocks.clear();
333 :
334 : int i = 0;
335 322 : for (const auto & b : getMooseMesh().meshSubdomains())
336 210 : _blocks[b] = i++;
337 :
338 112 : _n_block_bins = _blocks.size();
339 112 : }
340 :
341 : void
342 111 : MoabSkinner::createMOABElems()
343 : {
344 : // Clear prior results
345 : _id_to_elem_handles.clear();
346 :
347 : std::map<dof_id_type, moab::EntityHandle> node_id_to_handle;
348 :
349 : double coords[3];
350 :
351 : // Save all the node information
352 25832 : for (const auto & node : getMooseMesh().getMesh().node_ptr_range())
353 : {
354 : // Fetch coords (and scale to correct units)
355 25610 : coords[0] = _scaling * (*node)(0);
356 25610 : coords[1] = _scaling * (*node)(1);
357 25610 : coords[2] = _scaling * (*node)(2);
358 :
359 : // Add node to MOAB database and get handle
360 25610 : moab::EntityHandle ent(0);
361 25610 : check(_moab->create_vertex(coords, ent));
362 :
363 : // Save mapping of libMesh IDs to MOAB vertex handles
364 25610 : node_id_to_handle[node->id()] = ent;
365 111 : }
366 :
367 : moab::Range all_elems;
368 :
369 : // Iterate over elements in the mesh
370 135566 : for (const auto & elem : getMooseMesh().getMesh().active_element_ptr_range())
371 : {
372 67673 : auto nodeSets = getTetSets(elem->type());
373 :
374 : // Get the connectivity
375 : std::vector<dof_id_type> conn_libmesh;
376 67672 : elem->connectivity(0, libMesh::IOPackage::VTK, conn_libmesh);
377 :
378 : // Loop over sub tets
379 157968 : for (const auto & nodeSet : nodeSets)
380 : {
381 : // Set MOAB connectivity
382 90296 : std::vector<moab::EntityHandle> conn(NODES_PER_MOAB_TET);
383 451480 : for (unsigned int i = 0; i < NODES_PER_MOAB_TET; ++i)
384 : {
385 : // Get the elem node index of the ith node of the sub-tet
386 361184 : unsigned int nodeIndex = nodeSet.at(i);
387 361184 : conn[i] = node_id_to_handle[conn_libmesh.at(nodeIndex)];
388 : }
389 :
390 : // Create an element in MOAB database
391 90296 : moab::EntityHandle ent(0);
392 90296 : check(_moab->create_element(moab::MBTET, conn.data(), NODES_PER_MOAB_TET, ent));
393 :
394 : // Save mapping between libMesh ids and moab handles
395 90296 : auto id = elem->id();
396 90296 : if (_id_to_elem_handles.find(id) == _id_to_elem_handles.end())
397 135344 : _id_to_elem_handles[id] = std::vector<moab::EntityHandle>();
398 :
399 90296 : _id_to_elem_handles[id].push_back(ent);
400 :
401 : // Save the handle for adding to entity sets
402 90296 : all_elems.insert(ent);
403 90296 : }
404 67782 : }
405 :
406 : // Add the elems to the full meshset
407 110 : check(_moab->add_entities(_all_tets, all_elems));
408 :
409 : // Save the first elem
410 110 : offset = all_elems.front();
411 110 : }
412 :
413 : const std::vector<std::vector<unsigned int>> &
414 67673 : MoabSkinner::getTetSets(ElemType type) const
415 : {
416 67673 : if (type == TET4)
417 64440 : return _tet4_nodes;
418 3233 : else if (type == TET10)
419 3232 : return _tet10_nodes;
420 : else
421 1 : mooseError("The MoabSkinner can only be used with a tetrahedral [Mesh]!");
422 : }
423 :
424 : void
425 111 : MoabSkinner::createTags()
426 : {
427 : // First some built-in MOAB tag types
428 111 : check(_moab->tag_get_handle(GEOM_DIMENSION_TAG_NAME,
429 : 1,
430 : moab::MB_TYPE_INTEGER,
431 111 : geometry_dimension_tag,
432 : moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
433 :
434 111 : check(_moab->tag_get_handle(GLOBAL_ID_TAG_NAME,
435 : 1,
436 : moab::MB_TYPE_INTEGER,
437 111 : id_tag,
438 : moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
439 :
440 111 : check(_moab->tag_get_handle(CATEGORY_TAG_NAME,
441 : CATEGORY_TAG_SIZE,
442 : moab::MB_TYPE_OPAQUE,
443 111 : category_tag,
444 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
445 :
446 111 : check(_moab->tag_get_handle(NAME_TAG_NAME,
447 : NAME_TAG_SIZE,
448 : moab::MB_TYPE_OPAQUE,
449 111 : name_tag,
450 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
451 :
452 : // Some tags needed for DagMC
453 111 : check(_moab->tag_get_handle("FACETING_TOL",
454 : 1,
455 : moab::MB_TYPE_DOUBLE,
456 111 : faceting_tol_tag,
457 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
458 :
459 111 : check(_moab->tag_get_handle("GEOMETRY_RESABS",
460 : 1,
461 : moab::MB_TYPE_DOUBLE,
462 111 : geometry_resabs_tag,
463 : moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT));
464 :
465 : // Set the values for DagMC faceting / geometry tolerance tags on the mesh entity set
466 111 : check(_moab->tag_set_data(faceting_tol_tag, &_all_tets, 1, &_faceting_tol));
467 111 : check(_moab->tag_set_data(geometry_resabs_tag, &_all_tets, 1, &_geom_tol));
468 111 : }
469 :
470 : void
471 632 : MoabSkinner::createGroup(const unsigned int & id,
472 : const std::string & name,
473 : moab::EntityHandle & group_set)
474 : {
475 632 : check(_moab->create_meshset(moab::MESHSET_SET, group_set));
476 1264 : setTags(group_set, name, "Group", id, 4);
477 632 : }
478 :
479 : void
480 1640 : MoabSkinner::createVol(const unsigned int & id,
481 : moab::EntityHandle & volume_set,
482 : moab::EntityHandle group_set)
483 : {
484 1640 : check(_moab->create_meshset(moab::MESHSET_SET, volume_set));
485 :
486 3280 : setTags(volume_set, "", "Volume", id, 3);
487 :
488 : // Add the volume to group
489 1640 : check(_moab->add_entities(group_set, &volume_set, 1));
490 1640 : }
491 :
492 : void
493 3796 : MoabSkinner::createSurf(const unsigned int & id,
494 : moab::EntityHandle & surface_set,
495 : moab::Range & faces,
496 : const std::vector<VolData> & voldata)
497 : {
498 : // Create meshset
499 3796 : check(_moab->create_meshset(moab::MESHSET_SET, surface_set));
500 :
501 : // Set tags
502 7592 : setTags(surface_set, "", "Surface", id, 2);
503 :
504 : // Add tris to the surface
505 3796 : check(_moab->add_entities(surface_set, faces));
506 :
507 : // Create entry in map
508 7592 : surfsToVols[surface_set] = std::vector<VolData>();
509 :
510 : // Add volume to list associated with this surface
511 9718 : for (const auto & data : voldata)
512 5922 : updateSurfData(surface_set, data);
513 3796 : }
514 :
515 : void
516 6779 : MoabSkinner::updateSurfData(moab::EntityHandle surface_set, const VolData & data)
517 : {
518 : // Add the surface to the volume set
519 6779 : check(_moab->add_parent_child(data.vol, surface_set));
520 :
521 : // Set the surfaces sense
522 6779 : gtt->set_sense(surface_set, data.vol, int(data.sense));
523 :
524 6779 : surfsToVols[surface_set].push_back(data);
525 6779 : }
526 :
527 : void
528 6068 : MoabSkinner::setTags(
529 : moab::EntityHandle ent, std::string name, std::string category, unsigned int id, int dim)
530 : {
531 : // Set the name tag
532 6068 : if (name != "")
533 1264 : setTagData(name_tag, ent, name, NAME_TAG_SIZE);
534 :
535 : // Set the category tag
536 6068 : if (category != "")
537 12136 : setTagData(category_tag, ent, category, CATEGORY_TAG_SIZE);
538 :
539 : // Set the dimension tag
540 6068 : setTagData(geometry_dimension_tag, ent, &dim);
541 :
542 : // Set the id tag
543 6068 : setTagData(id_tag, ent, &id);
544 6068 : }
545 :
546 : void
547 6700 : MoabSkinner::setTagData(moab::Tag tag, moab::EntityHandle ent, std::string data, unsigned int SIZE)
548 : {
549 6700 : auto namebuf = new char[SIZE];
550 6700 : memset(namebuf, '\0', SIZE); // fill C char array with null
551 6700 : strncpy(namebuf, data.c_str(), SIZE - 1);
552 6700 : check(_moab->tag_set_data(tag, &ent, 1, namebuf));
553 6700 : delete[] namebuf;
554 6700 : }
555 :
556 : void
557 12136 : MoabSkinner::setTagData(moab::Tag tag, moab::EntityHandle ent, void * data)
558 : {
559 12136 : check(_moab->tag_set_data(tag, &ent, 1, data));
560 12136 : }
561 :
562 : unsigned int
563 78 : MoabSkinner::nBins() const
564 : {
565 78 : return _n_block_bins * _n_density_bins * _n_temperature_bins;
566 : }
567 :
568 : void
569 46 : MoabSkinner::sortElemsByResults()
570 : {
571 46 : _elem_bins.clear();
572 46 : _elem_bins.resize(nBins());
573 :
574 : // accumulate information for printing diagnostics
575 46 : std::vector<unsigned int> n_block_hits(_n_block_bins, 0);
576 46 : std::vector<unsigned int> n_temp_hits(_n_temperature_bins, 0);
577 46 : std::vector<unsigned int> n_density_hits(_n_density_bins, 0);
578 :
579 25838 : for (unsigned int e = 0; e < getMooseMesh().nElem(); ++e)
580 : {
581 25796 : const Elem * const elem = getMooseMesh().queryElemPtr(e);
582 25796 : if (!elem)
583 0 : continue;
584 :
585 : // bin by subdomain ID
586 25796 : auto iMat = getSubdomainBin(elem);
587 25796 : n_block_hits[iMat] += 1;
588 :
589 : // bin by density
590 25796 : auto iDenBin = getDensityBin(elem);
591 25794 : n_density_hits[iDenBin] += 1;
592 :
593 : // bin by temperature
594 25794 : auto iBin = getTemperatureBin(elem);
595 25792 : n_temp_hits[iBin] += 1;
596 :
597 : // Sort elem into a bin
598 25792 : auto iSortBin = getBin(iBin, iDenBin, iMat);
599 25792 : _elem_bins.at(iSortBin).insert(elem->id());
600 : }
601 :
602 42 : if (_verbose)
603 : {
604 40 : VariadicTable<unsigned int, std::string, unsigned int> vtt({"Bin", "Range (K)", "# Elems"});
605 40 : VariadicTable<unsigned int, std::string, unsigned int> vtd({"Bin", "Range (kg/m3)", "# Elems"});
606 :
607 186 : for (unsigned int i = 0; i < _n_temperature_bins; ++i)
608 146 : vtt.addRow(i,
609 438 : std::to_string(_temperature_bin_bounds[i]) + " to " +
610 292 : std::to_string(_temperature_bin_bounds[i + 1]),
611 146 : n_temp_hits[i]);
612 :
613 130 : for (unsigned int i = 0; i < _n_density_bins; ++i)
614 90 : vtd.addRow(i,
615 270 : std::to_string(_density_bin_bounds[i]) + " to " +
616 180 : std::to_string(_density_bin_bounds[i + 1]),
617 90 : n_density_hits[i]);
618 :
619 40 : _console << "\nMapping of Elements to Temperature Bins:" << std::endl;
620 40 : vtt.print(_console);
621 40 : _console << std::endl;
622 :
623 40 : if (_bin_by_density)
624 : {
625 13 : _console << "\n\nMapping of Elements to Density Bins:" << std::endl;
626 13 : vtd.print(_console);
627 13 : _console << std::endl;
628 : }
629 40 : }
630 42 : }
631 :
632 : unsigned int
633 51650 : MoabSkinner::getTemperatureBin(const Elem * const elem) const
634 : {
635 51650 : auto dof = elem->dof_number(_fe_problem.getAuxiliarySystem().number(), _temperature_var_num, 0);
636 51650 : auto value = (*_serialized_solution)(dof);
637 :
638 : // TODO: add option to truncate instead
639 51650 : if ((_temperature_min - value) > BIN_TOLERANCE)
640 1 : mooseError("Variable '",
641 1 : _temperature_name,
642 : "' has value below minimum range of bins. "
643 : "Please decrease 'temperature_min'.\n\n"
644 : " value: ",
645 : value,
646 : "\n temperature_min: ",
647 : _temperature_min);
648 :
649 51649 : if ((value - _temperature_max) > BIN_TOLERANCE)
650 1 : mooseError("Variable '",
651 1 : _temperature_name,
652 : "' has value above maximum range of bins. "
653 : "Please increase 'temperature_max'.\n\n"
654 : " value: ",
655 : value,
656 : "\n temperature_max: ",
657 : _temperature_max);
658 :
659 51648 : return bin_utility::linearBin(value, _temperature_bin_bounds);
660 : }
661 :
662 : unsigned int
663 50036 : MoabSkinner::getDensityBin(const Elem * const elem) const
664 : {
665 50036 : if (!_bin_by_density)
666 : return 0;
667 :
668 30588 : auto dof = elem->dof_number(_fe_problem.getAuxiliarySystem().number(), _density_var_num, 0);
669 30588 : auto value = (*_serialized_solution)(dof);
670 :
671 : // TODO: add option to truncate instead
672 30588 : if ((_density_min - value) > BIN_TOLERANCE)
673 1 : mooseError("Variable '",
674 1 : _density_name,
675 : "' has value below minimum range of bins. "
676 : "Please decrease 'density_min'.\n\n"
677 : " value: ",
678 : value,
679 : "\n density_min: ",
680 1 : _density_min);
681 :
682 30587 : if ((value - _density_max) > BIN_TOLERANCE)
683 1 : mooseError("Variable '",
684 1 : _density_name,
685 : "' has value above maximum range of bins. "
686 : "Please increase 'density_max'.\n\n"
687 : " value: ",
688 : value,
689 : "\n density_max: ",
690 1 : _density_max);
691 :
692 30586 : return bin_utility::linearBin(value, _density_bin_bounds);
693 : }
694 :
695 : std::string
696 600 : MoabSkinner::materialName(const unsigned int & block,
697 : const unsigned int & density,
698 : const unsigned int & temp) const
699 : {
700 600 : return "mat:" + _material_names.at(block);
701 : }
702 :
703 : void
704 42 : MoabSkinner::findSurfaces()
705 : {
706 : // Find all neighbours in mesh
707 42 : getMooseMesh().getMesh().find_neighbors();
708 :
709 : // Counter for volumes
710 42 : unsigned int vol_id = 0;
711 :
712 : // Counter for surfaces
713 42 : unsigned int surf_id = 0;
714 :
715 : // Loop over material bins
716 124 : for (unsigned int iMat = 0; iMat < _n_block_bins; iMat++)
717 : {
718 : // Loop over density bins
719 246 : for (unsigned int iDen = 0; iDen < _n_density_bins; iDen++)
720 : {
721 : // Loop over temperature bins
722 764 : for (unsigned int iVar = 0; iVar < _n_temperature_bins; iVar++)
723 : {
724 : // Update material name
725 600 : auto updated_mat_name = materialName(iMat, iDen, iVar);
726 :
727 : // Create a material group
728 600 : int iSortBin = getBin(iVar, iDen, iMat);
729 :
730 : // For DagMC to fill a cell with a material, we first create a group
731 : // with that name, and then assign it with createVol (called inside findSurface)
732 : moab::EntityHandle group_set;
733 600 : unsigned int group_id = iSortBin + 1;
734 600 : createGroup(group_id, updated_mat_name, group_set);
735 :
736 : // Sort elems in this mat-density-temp bin into local regions
737 : std::vector<moab::Range> regions;
738 1800 : groupLocalElems(_elem_bins.at(iSortBin), regions);
739 :
740 : // Loop over all regions and find surfaces
741 2210 : for (const auto & region : regions)
742 : {
743 : moab::EntityHandle volume_set;
744 1610 : findSurface(region, group_set, vol_id, surf_id, volume_set);
745 : }
746 600 : }
747 : }
748 : }
749 :
750 42 : if (_build_graveyard)
751 30 : buildGraveyard(vol_id, surf_id);
752 :
753 42 : if (_set_implicit_complement_material)
754 : {
755 : moab::EntityHandle comp_group;
756 2 : unsigned int comp_id = nBins() + 1 + _build_graveyard;
757 2 : createGroup(comp_id, _implicit_complement_group_name, comp_group);
758 2 : moab::EntityHandle arbitray_volume = 0;
759 2 : for (const auto & surf_pair : surfsToVols)
760 : {
761 : const auto & vols = surf_pair.second;
762 2 : arbitray_volume = vols.front().vol;
763 2 : break;
764 : }
765 2 : check(_moab->add_entities(comp_group, &arbitray_volume, 1));
766 : }
767 :
768 : // Write MOAB volume and/or skin meshes to file
769 42 : write();
770 42 : }
771 :
772 : void
773 42 : MoabSkinner::write()
774 : {
775 : // Only write to file on root process
776 42 : if (processor_id() != 0)
777 0 : return;
778 :
779 42 : std::string extension = std::to_string(_n_write) + ".h5m";
780 :
781 42 : if (_output_skins)
782 : {
783 : // Generate list of surfaces to write
784 : std::vector<moab::EntityHandle> surfs;
785 3271 : for (const auto & itsurf : surfsToVols)
786 3260 : surfs.push_back(itsurf.first);
787 :
788 11 : std::string filename = "moab_skins_" + extension;
789 :
790 11 : if (_verbose)
791 11 : _console << "Writing MOAB skins to " << filename << "...";
792 :
793 11 : check(_moab->write_mesh(filename.c_str(), surfs.data(), surfs.size()));
794 11 : }
795 :
796 42 : if (_output_full)
797 : {
798 1 : std::string filename = "moab_mesh_" + extension;
799 :
800 1 : if (_verbose)
801 1 : _console << "Writing MOAB mesh to " << filename << std::endl;
802 :
803 1 : check(_moab->write_mesh(filename.c_str()));
804 : }
805 :
806 42 : _n_write++;
807 : }
808 :
809 : void
810 600 : MoabSkinner::groupLocalElems(std::set<dof_id_type> elems, std::vector<moab::Range> & localElems)
811 : {
812 2210 : while (!elems.empty())
813 : {
814 :
815 : // Create a new local range of moab handles
816 : moab::Range local;
817 :
818 : // Retrieve and remove the fisrt elem
819 : auto it = elems.begin();
820 1610 : dof_id_type next = *it;
821 1610 : elems.erase(it);
822 :
823 : std::set<dof_id_type> neighbors;
824 1610 : neighbors.insert(next);
825 :
826 6378 : while (!neighbors.empty())
827 : {
828 :
829 : std::set<dof_id_type> new_neighbors;
830 :
831 : // Loop over all the new neighbors
832 30560 : for (auto & next : neighbors)
833 : {
834 :
835 : // Get the MOAB handles, and add to local set
836 : // (May be more than one if this libMesh elem has sub-tetrahedra)
837 25792 : if (_id_to_elem_handles.find(next) == _id_to_elem_handles.end())
838 0 : mooseError("No entity handles found for libmesh id.");
839 :
840 25792 : std::vector<moab::EntityHandle> ents = _id_to_elem_handles[next];
841 62896 : for (const auto ent : ents)
842 37104 : local.insert(ent);
843 :
844 : // Get the libMesh element
845 25792 : Elem & elem = getMooseMesh().getMesh().elem_ref(next);
846 :
847 : // How many nearest neighbors (general element)?
848 : unsigned int NN = elem.n_neighbors();
849 :
850 : // Loop over neighbors
851 128960 : for (unsigned int i = 0; i < NN; i++)
852 : {
853 : const Elem * nnptr = elem.neighbor_ptr(i);
854 : // If on boundary, some may be null ptrs
855 103168 : if (nnptr == nullptr)
856 11744 : continue;
857 :
858 91424 : dof_id_type idnn = nnptr->id();
859 :
860 : // Select only those that are in the current bin
861 91424 : if (elems.find(idnn) != elems.end())
862 : {
863 24182 : new_neighbors.insert(idnn);
864 : // Remove from those still available
865 : elems.erase(idnn);
866 : }
867 : }
868 25792 : }
869 :
870 : // Found all the new neighbors, done with current set.
871 : neighbors = new_neighbors;
872 : }
873 :
874 : // Save this moab range of local neighbors
875 1610 : localElems.push_back(local);
876 : }
877 600 : }
878 :
879 : void
880 46 : MoabSkinner::reset()
881 : {
882 46 : _moab.reset(new moab::Core());
883 46 : skinner.reset(new moab::Skinner(_moab.get()));
884 46 : gtt.reset(new moab::GeomTopoTool(_moab.get()));
885 :
886 : // Clear entity set maps
887 : surfsToVols.clear();
888 46 : }
889 :
890 : unsigned int
891 40128 : MoabSkinner::getBin(const unsigned int & i_temp,
892 : const unsigned int & i_density,
893 : const unsigned int & i_block) const
894 : {
895 40128 : return _n_temperature_bins * (_n_density_bins * i_block + i_density) + i_temp;
896 : }
897 :
898 : void
899 1610 : MoabSkinner::findSurface(const moab::Range & region,
900 : moab::EntityHandle group,
901 : unsigned int & vol_id,
902 : unsigned int & surf_id,
903 : moab::EntityHandle & volume_set)
904 : {
905 : // Create a volume set
906 1610 : vol_id++;
907 1610 : createVol(vol_id, volume_set, group);
908 :
909 : // Find surfaces from these regions
910 : moab::Range tris; // The tris of the surfaces
911 : moab::Range rtris; // The tris which are reversed with respect to their surfaces
912 1610 : skinner->find_skin(0, region, false, tris, &rtris);
913 :
914 : // Create surface sets for the forwards tris
915 1610 : VolData vdata = {volume_set, Sense::FORWARDS};
916 1610 : createSurfaces(tris, vdata, surf_id);
917 :
918 : // Create surface sets for the reversed tris
919 1610 : vdata.sense = Sense::BACKWARDS;
920 1610 : createSurfaces(rtris, vdata, surf_id);
921 1610 : }
922 :
923 : void
924 3220 : MoabSkinner::createSurfaces(moab::Range & faces, VolData & voldata, unsigned int & surf_id)
925 : {
926 3220 : if (faces.empty())
927 : return;
928 :
929 : // Loop over the surfaces we have already created
930 1793957 : for (const auto & surfpair : surfsToVols)
931 : {
932 : // Local copies of surf/vols
933 1792347 : moab::EntityHandle surf = surfpair.first;
934 1792347 : std::vector<VolData> vols = surfpair.second;
935 :
936 : // First get the entities in this surface
937 : moab::Range tris;
938 1792347 : check(_moab->get_entities_by_handle(surf, tris));
939 :
940 : // Find any tris that live in both surfs
941 1792347 : moab::Range overlap = moab::intersect(tris, faces);
942 1792347 : if (!overlap.empty())
943 : {
944 : // Check if the tris are a subset or the entire surf
945 2983 : if (tris.size() == overlap.size())
946 : {
947 : // Whole surface -> Just need to update the volume relationships
948 857 : updateSurfData(surf, voldata);
949 : }
950 : else
951 : {
952 : // If overlap is subset, subtract shared tris from this surface and create a new shared
953 : // surface
954 2126 : check(_moab->remove_entities(surf, overlap));
955 :
956 : // Append our new volume to the list that share this surf
957 2126 : vols.push_back(voldata);
958 :
959 : // Create a new shared surface
960 : moab::EntityHandle shared_surf;
961 2126 : surf_id++;
962 2126 : createSurf(surf_id, shared_surf, overlap, vols);
963 : }
964 :
965 : // Subtract from the input list
966 10993 : for (auto & shared : overlap)
967 8010 : faces.erase(shared);
968 :
969 2983 : if (faces.empty())
970 : break;
971 : }
972 1792347 : }
973 :
974 3138 : if (!faces.empty())
975 : {
976 : moab::EntityHandle surface_set;
977 1610 : std::vector<VolData> voldatavec(1, voldata);
978 1610 : surf_id++;
979 1610 : createSurf(surf_id, surface_set, faces, voldatavec);
980 1610 : }
981 : }
982 :
983 : void
984 30 : MoabSkinner::buildGraveyard(unsigned int & vol_id, unsigned int & surf_id)
985 : {
986 : // Create the graveyard group
987 : moab::EntityHandle graveyard;
988 30 : unsigned int id = nBins() + 1;
989 30 : createGroup(id, "mat:Graveyard", graveyard);
990 :
991 : // Create a volume set
992 : moab::EntityHandle volume_set;
993 30 : createVol(++vol_id, volume_set, graveyard);
994 :
995 : // Set up for the volume data to pass to surfs
996 30 : VolData vdata = {volume_set, Sense::FORWARDS};
997 :
998 : // Find a bounding box
999 30 : BoundingBox bbox = MeshTools::create_bounding_box(getMooseMesh().getMesh());
1000 :
1001 : // Build the two cubic surfaces defining the graveyard
1002 30 : createSurfaceFromBox(
1003 : bbox, vdata, surf_id, false /* normals point into box */, _graveyard_scale_inner);
1004 30 : createSurfaceFromBox(
1005 : bbox, vdata, surf_id, true /* normals point out of box */, _graveyard_scale_outer);
1006 30 : }
1007 :
1008 : void
1009 60 : MoabSkinner::createSurfaceFromBox(const BoundingBox & box,
1010 : const VolData & voldata,
1011 : unsigned int & surf_id,
1012 : bool normalout,
1013 : const Real & factor)
1014 : {
1015 60 : std::vector<moab::EntityHandle> vert_handles = createNodesFromBox(box, factor);
1016 :
1017 : // Create the tris in 4 groups of 3 (4 open tetrahedra)
1018 : moab::Range tris;
1019 60 : createCornerTris(vert_handles, 0, 1, 2, 4, normalout, tris);
1020 60 : createCornerTris(vert_handles, 3, 2, 1, 7, normalout, tris);
1021 60 : createCornerTris(vert_handles, 6, 4, 2, 7, normalout, tris);
1022 60 : createCornerTris(vert_handles, 5, 1, 4, 7, normalout, tris);
1023 :
1024 : moab::EntityHandle surface_set;
1025 60 : std::vector<VolData> voldatavec(1, voldata);
1026 60 : surf_id++;
1027 60 : createSurf(surf_id, surface_set, tris, voldatavec);
1028 120 : }
1029 :
1030 : std::vector<moab::EntityHandle>
1031 60 : MoabSkinner::createNodesFromBox(const BoundingBox & box, const Real & factor) const
1032 : {
1033 : std::vector<moab::EntityHandle> vert_handles;
1034 :
1035 : // Fetch the vertices of the box
1036 60 : auto verts = geom_utils::boxCorners(box, factor);
1037 :
1038 : // Array to represent a coord in MOAB
1039 : double coord[3];
1040 :
1041 : // Create the vertices in MOAB and get the handles
1042 540 : for (const auto & vert : verts)
1043 : {
1044 480 : coord[0] = vert(0) * _scaling;
1045 480 : coord[1] = vert(1) * _scaling;
1046 480 : coord[2] = vert(2) * _scaling;
1047 :
1048 : moab::EntityHandle ent;
1049 480 : check(_moab->create_vertex(coord, ent));
1050 480 : vert_handles.push_back(ent);
1051 : }
1052 :
1053 60 : return vert_handles;
1054 60 : }
1055 :
1056 : void
1057 240 : MoabSkinner::createCornerTris(const std::vector<moab::EntityHandle> & verts,
1058 : unsigned int corner,
1059 : unsigned int v1,
1060 : unsigned int v2,
1061 : unsigned int v3,
1062 : bool normalout,
1063 : moab::Range & surface_tris)
1064 : {
1065 : // Create 3 tris stemming from one corner (i.e. an open tetrahedron)
1066 : // Assume first is the central corner, and the others are labelled clockwise looking down on the
1067 : // corner
1068 240 : unsigned int indices[3] = {v1, v2, v3};
1069 :
1070 : // Create each tri by a cyclic permutation of indices
1071 : // Values of i1, i2 in the loop: 0,1; 1,2; 2;0
1072 960 : for (unsigned int i = 0; i < 3; i++)
1073 : {
1074 720 : int i1 = indices[i % 3];
1075 720 : int i2 = indices[(i + 1) % 3];
1076 720 : if (normalout) // anti-clockwise: normal points outwards
1077 360 : surface_tris.insert(createTri(verts, corner, i2, i1));
1078 : else // clockwise: normal points inwards
1079 360 : surface_tris.insert(createTri(verts, corner, i1, i2));
1080 : }
1081 240 : }
1082 :
1083 : moab::EntityHandle
1084 720 : MoabSkinner::createTri(const std::vector<moab::EntityHandle> & vertices,
1085 : unsigned int v1,
1086 : unsigned int v2,
1087 : unsigned int v3)
1088 : {
1089 : moab::EntityHandle triangle;
1090 720 : moab::EntityHandle connectivity[3] = {vertices[v1], vertices[v2], vertices[v3]};
1091 720 : check(_moab->create_element(moab::MBTRI, connectivity, 3, triangle));
1092 720 : return triangle;
1093 : }
1094 :
1095 : void
1096 18 : MoabSkinner::setGraveyard(bool build)
1097 : {
1098 18 : if (build != _build_graveyard)
1099 : {
1100 2 : std::string original = _build_graveyard ? "true" : "false";
1101 2 : std::string change = _build_graveyard ? "false" : "true";
1102 1 : mooseWarning("Overriding graveyard setting from ",
1103 : original,
1104 : " to ",
1105 : change,
1106 : ".\n"
1107 : "To hide this warning, set 'build_graveyard = ",
1108 : change,
1109 : "'");
1110 : }
1111 :
1112 17 : _build_graveyard = build;
1113 17 : }
1114 : #endif
|
