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