Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CoarsenBlockGenerator.h"
11 : #include "MooseMeshUtils.h"
12 : #include "MeshCoarseningUtils.h"
13 : #include "MeshBaseDiagnosticsUtils.h"
14 :
15 : #include "libmesh/elem.h"
16 : #include "libmesh/mesh_modification.h"
17 : #include "CastUniquePointer.h"
18 :
19 : registerMooseObject("MooseApp", CoarsenBlockGenerator);
20 :
21 : InputParameters
22 3468 : CoarsenBlockGenerator::validParams()
23 : {
24 3468 : InputParameters params = MeshGenerator::validParams();
25 :
26 6936 : params.addClassDescription("Mesh generator which coarsens one or more blocks in an existing "
27 : "mesh. The coarsening algorithm works best for regular meshes.");
28 13872 : params.addRequiredParam<MeshGeneratorName>("input", "Input mesh to coarsen");
29 13872 : params.addRequiredParam<std::vector<SubdomainName>>("block",
30 : "The list of blocks to be coarsened");
31 13872 : params.addRequiredParam<std::vector<unsigned int>>(
32 : "coarsening",
33 : "Maximum amount of times to coarsen elements in each block. See 'block' for indexing");
34 13872 : params.addRequiredParam<Point>("starting_point",
35 : "A point inside the element to start the coarsening from");
36 :
37 : // This is a heuristic to be able to coarsen inside blocks that are not uniformly refined
38 20808 : params.addRangeCheckedParam<Real>(
39 : "maximum_volume_ratio",
40 : 2,
41 : "maximum_volume_ratio > 0",
42 : "Maximum allowed volume ratio between two fine elements to propagate "
43 : "the coarsening front through a side");
44 10404 : params.addParam<bool>(
45 : "verbose",
46 6936 : false,
47 : "Whether to make the mesh generator output details of its actions on the console");
48 6936 : params.addParam<bool>("check_for_non_conformal_output_mesh",
49 6936 : true,
50 : "Whether to check the entire mesh for non-conformal nodes indicating that "
51 : "the coarsening operation has failed to produce a conformal mesh");
52 3468 : return params;
53 0 : }
54 :
55 205 : CoarsenBlockGenerator::CoarsenBlockGenerator(const InputParameters & parameters)
56 : : MeshGenerator(parameters),
57 205 : _input(getMesh("input")),
58 410 : _block(getParam<std::vector<SubdomainName>>("block")),
59 410 : _coarsening(getParam<std::vector<unsigned int>>("coarsening")),
60 410 : _starting_point(getParam<Point>("starting_point")),
61 410 : _max_vol_ratio(getParam<Real>("maximum_volume_ratio")),
62 410 : _verbose(getParam<bool>("verbose")),
63 615 : _check_output_mesh_for_nonconformality(getParam<bool>("check_for_non_conformal_output_mesh"))
64 : {
65 205 : if (_block.size() != _coarsening.size())
66 6 : paramError("coarsening", "The blocks and coarsening parameter vectors should be the same size");
67 202 : }
68 :
69 : std::unique_ptr<MeshBase>
70 102 : CoarsenBlockGenerator::generate()
71 : {
72 : // Get the list of block ids from the block names
73 : const auto block_ids =
74 204 : MooseMeshUtils::getSubdomainIDs(*_input, getParam<std::vector<SubdomainName>>("block"));
75 102 : const std::set<SubdomainID> block_ids_set(block_ids.begin(), block_ids.end());
76 :
77 : // Check that the block ids/names exist in the mesh
78 102 : if (!_input->preparation().has_cached_elem_data)
79 19 : _input->cache_elem_data();
80 102 : std::set<SubdomainID> mesh_blocks;
81 102 : _input->subdomain_ids(mesh_blocks);
82 :
83 367 : for (std::size_t i = 0; i < block_ids.size(); ++i)
84 268 : if (!MooseMeshUtils::hasSubdomainID(*_input, block_ids[i]))
85 6 : paramError("block",
86 : "The block '",
87 6 : getParam<std::vector<SubdomainName>>("block")[i],
88 : "' was not found within the mesh");
89 :
90 : // Error if it has not been implemented for this element type
91 29500 : for (const auto & elem : _input->active_subdomain_set_elements_ptr_range(block_ids_set))
92 : // Only types implemented
93 29401 : if (elem->type() != QUAD4 && elem->type() != HEX8)
94 0 : paramError("block",
95 0 : "The input mesh contains an unsupported element type '" +
96 0 : Moose::stringify(elem->type()) + "' for coarsening in block " +
97 99 : std::to_string(elem->subdomain_id()));
98 :
99 : // Take ownership of the mesh
100 99 : std::unique_ptr<MeshBase> mesh = std::move(_input);
101 99 : if (!mesh->is_serial())
102 0 : paramError("input", "Input mesh must not be distributed");
103 :
104 : // Find the element to start from
105 99 : auto start_elem = (*mesh->sub_point_locator())(_starting_point);
106 :
107 : // Check that the starting element choice: in the block with the most coarsening requested
108 99 : if (!start_elem)
109 0 : paramError("starting_point", "No element was found at that point");
110 99 : unsigned int max_c = *std::max_element(_coarsening.begin(), _coarsening.end());
111 358 : for (const auto i : index_range(block_ids))
112 262 : if (block_ids[i] == start_elem->subdomain_id() && _coarsening[i] != max_c)
113 6 : mooseError("The starting element must be in the block set to be coarsened the most.\n"
114 : "Starting element is in block ",
115 3 : start_elem->subdomain_id(),
116 : " set to be coarsened ",
117 3 : _coarsening[i],
118 : " times but the max coarsening required is ",
119 : max_c);
120 :
121 : // Determine how many times the coarsening will be used
122 96 : if (max_c > 0 && !mesh->is_prepared())
123 : // we prepare for use to make sure the neighbors have been found
124 16 : mesh->prepare_for_use();
125 :
126 96 : auto mesh_ptr = recursiveCoarsen(block_ids, mesh, _coarsening, max_c, /*step=*/0);
127 :
128 : // element neighbors are not valid
129 96 : if (max_c > 0)
130 96 : mesh_ptr->unset_is_prepared();
131 :
132 : // flip elements as we were not careful to build them with a positive volume
133 96 : MeshTools::Modification::orient_elements(*mesh_ptr);
134 :
135 : // check that we are not returning a non-conformal mesh
136 96 : if (_check_output_mesh_for_nonconformality)
137 : {
138 96 : mesh_ptr->prepare_for_use();
139 96 : unsigned int num_nonconformal_nodes = 0;
140 96 : MeshBaseDiagnosticsUtils::checkNonConformalMesh(
141 96 : mesh_ptr, _console, 10, TOLERANCE, num_nonconformal_nodes);
142 96 : if (num_nonconformal_nodes)
143 0 : mooseError("Coarsened mesh has non-conformal nodes. The coarsening process likely failed to "
144 0 : "form a uniform paving of coarsened elements. Number of non-conformal nodes: " +
145 0 : Moose::stringify(num_nonconformal_nodes));
146 : }
147 192 : return mesh_ptr;
148 96 : }
149 :
150 : std::unique_ptr<MeshBase>
151 208 : CoarsenBlockGenerator::recursiveCoarsen(const std::vector<subdomain_id_type> & block_ids,
152 : std::unique_ptr<MeshBase> & mesh,
153 : const std::vector<unsigned int> & coarsening,
154 : const unsigned int max,
155 : unsigned int coarse_step)
156 : {
157 208 : if (coarse_step == max)
158 96 : return dynamic_pointer_cast<MeshBase>(mesh);
159 :
160 : // Elements should know their neighbors
161 112 : if (!mesh->is_prepared())
162 0 : mesh->prepare_for_use();
163 :
164 : // We wont be modifying the starting mesh for simplicity, we will make a copy and return that
165 112 : std::unique_ptr<MeshBase> mesh_return;
166 112 : int max_num_coarsened = -1;
167 :
168 112 : const auto base_start_elem = (*mesh->sub_point_locator())(_starting_point);
169 :
170 : // Try every node as the 'center' point of a coarsened element
171 784 : for (const auto & start_node_index : base_start_elem->node_index_range())
172 : {
173 672 : if (_verbose)
174 576 : _console << "Step " << coarse_step + 1 << " coarsening attempt #" << start_node_index
175 576 : << "\nUsing node " << *base_start_elem->node_ptr(start_node_index)
176 576 : << " as the interior node of the coarse element." << std::endl;
177 :
178 : // Make a copy of the mesh in case the initial node choice was bad
179 672 : auto mesh_copy = mesh->clone();
180 :
181 : // We will only have a single starting element for now. If there are non-connected components,
182 : // we will need to have a starting element-node pair in every component.
183 672 : auto start_elem = mesh_copy->elem_ptr(base_start_elem->id());
184 : mooseAssert(start_elem, "Should have a real elem pointer");
185 : mooseAssert(start_elem->active(), "Starting element must be active");
186 :
187 672 : auto start_node = start_elem->node_ptr(start_node_index);
188 : mooseAssert(start_node, "Starting node should exist");
189 :
190 : // Create comparator for ordering of candidates
191 93819 : auto cmp = [](std::pair<Elem *, Node *> a, std::pair<Elem *, Node *> b)
192 : {
193 : // Sweep direction
194 : // Potentially a user selectable parameter in the future
195 93819 : Point sorting_direction(1, 1, 1);
196 : const auto sorting =
197 93819 : (a.first->vertex_average() - b.first->vertex_average()) * sorting_direction;
198 93819 : if (MooseUtils::absoluteFuzzyGreaterThan(sorting, 0))
199 41083 : return true;
200 65482 : else if (MooseUtils::absoluteFuzzyEqual(sorting, 0) &&
201 65482 : MooseUtils::absoluteFuzzyGreaterThan((*a.second - *b.second) * sorting_direction, 0))
202 536 : return true;
203 : else
204 : // Sorting direction is orthogonal to the two pairs, rely on element ids
205 52200 : return a.first->id() > b.first->id();
206 : };
207 :
208 : // This set will keep track of all the 'fine elem' + 'coarse element interior node' pairs
209 : // we should attempt to form coarse element from
210 : // TODO: think about the implications of set vs vector. Set might grow to the entire mesh
211 : // due to sorting. Vector we could insert at the beginning and treat new candidates immediately
212 672 : std::set<std::pair<Elem *, Node *>, decltype(cmp)> candidate_pairs(cmp);
213 672 : candidate_pairs.insert(std::make_pair(start_elem, start_node));
214 :
215 : // Keep track of the coarse elements created
216 672 : std::set<Elem *> coarse_elems;
217 :
218 12040 : while (candidate_pairs.size() > 0)
219 : {
220 11368 : Elem * current_elem = candidate_pairs.begin()->first;
221 11368 : Node * interior_node = candidate_pairs.begin()->second;
222 : mooseAssert(current_elem, "Null candidate element pointer");
223 : mooseAssert(interior_node, "Null candidate node pointer");
224 11368 : const auto current_node_index = current_elem->get_node_index(interior_node);
225 : // just take any another node for now
226 : const auto ref_node =
227 11368 : current_elem->node_ptr(current_node_index == 0 ? 1 : current_node_index - 1);
228 : mooseAssert(ref_node, "Should have a real node pointer");
229 11368 : candidate_pairs.erase(candidate_pairs.begin());
230 :
231 11368 : const auto elem_type = current_elem->type();
232 :
233 : // Mid-edge nodes could be an option too for making coarse elements.
234 : // For a first implementation, we wont try to use them as near the edge we would need
235 : // a special treatment.
236 11368 : if (!current_elem->is_vertex(current_node_index))
237 4696 : continue;
238 :
239 : // We dont support coarsening libMesh h-refined meshes
240 11368 : if (current_elem->level() > 0)
241 0 : mooseError("H-refined meshes cannot be coarsened with this mesh generator. Use the "
242 : "[Adaptivity] block to coarsen them.");
243 :
244 : // Get the nodes to build a coarse element
245 11368 : std::vector<const Node *> tentative_coarse_nodes;
246 11368 : std::set<const Elem *> fine_elements_const;
247 11368 : bool success = MeshCoarseningUtils::getFineElementsFromInteriorNode(
248 : *interior_node, *ref_node, *current_elem, tentative_coarse_nodes, fine_elements_const);
249 :
250 : // For example, not enough fine elements around the node to build a coarse element
251 11368 : if (!success)
252 4056 : continue;
253 :
254 7312 : bool go_to_next_candidate = false;
255 : // If the fine elements are not all of the same type, we currently cannot coarsen
256 61872 : for (auto elem : fine_elements_const)
257 54560 : if (elem && elem->type() != elem_type)
258 0 : go_to_next_candidate = true;
259 :
260 : // We do not coarsen across subdomains for now
261 7312 : const auto common_subdomain_id = current_elem->subdomain_id();
262 61872 : for (auto elem : fine_elements_const)
263 54560 : if (elem && elem->subdomain_id() != common_subdomain_id)
264 2320 : go_to_next_candidate = true;
265 :
266 : // Check the coarse element nodes gathered
267 61872 : for (const auto & check_node : tentative_coarse_nodes)
268 54560 : if (check_node == nullptr)
269 0 : go_to_next_candidate = true;
270 :
271 7312 : if (go_to_next_candidate)
272 640 : continue;
273 :
274 : // We will likely delete the fine elements so we have to drop the const
275 255808 : auto cmp_elem = [](Elem * a, Elem * b) { return a->id() - b->id(); };
276 6672 : std::set<Elem *, decltype(cmp_elem)> fine_elements(cmp_elem);
277 56592 : for (const auto elem_ptr : fine_elements_const)
278 49920 : fine_elements.insert(mesh_copy->elem_ptr(elem_ptr->id()));
279 :
280 : // Form a parent, of a low order type as we only have the extreme vertex nodes
281 6672 : std::unique_ptr<Elem> parent = Elem::build(Elem::first_order_equivalent_type(elem_type));
282 6672 : parent->subdomain_id() = common_subdomain_id;
283 6672 : auto parent_ptr = mesh_copy->add_elem(parent.release());
284 6672 : coarse_elems.insert(parent_ptr);
285 :
286 : // Set the nodes to the coarse element
287 : // They were sorted previously in getFineElementFromInteriorNode
288 56592 : for (auto i : index_range(tentative_coarse_nodes))
289 49920 : parent_ptr->set_node(i, mesh_copy->node_ptr(tentative_coarse_nodes[i]->id()));
290 :
291 : // Gather targets / next candidates for the next element coarsening
292 : // Find the face neighbors, then look for the center node
293 44976 : for (const auto side_index : make_range(parent_ptr->n_sides()))
294 : {
295 : // Pick one of the coarse element nodes by that face
296 : // it should not matter which one, they are all vertex nodes of a fine element
297 : // that has a neighbor on the other side of the coarse element face
298 38304 : const auto coarse_node = parent_ptr->side_ptr(side_index)->node_ptr(0);
299 : mooseAssert(coarse_node,
300 : "We should have a node on coarse side " + std::to_string(side_index));
301 :
302 : // Find one of the fine elements next to the face, its neighbor on the other side
303 : // of the coarse face is the face neighbor we want
304 38304 : Elem * fine_el = nullptr;
305 139250 : for (const auto & fine_elem : fine_elements)
306 : {
307 139250 : bool found = false;
308 1094634 : for (const auto & fine_elem_node : fine_elem->node_ref_range())
309 993688 : if (MooseUtils::absoluteFuzzyEqual((*coarse_node - fine_elem_node).norm_sq(), 0))
310 : {
311 38304 : fine_el = fine_elem;
312 38304 : found = true;
313 38304 : break;
314 : }
315 139250 : if (found)
316 38304 : break;
317 : }
318 : mooseAssert(fine_el, "We should have found a fine element for the next candidate");
319 38304 : const Real fine_el_volume = fine_el->volume();
320 :
321 : // Get the element(s) on the other side of the coarse face
322 : // We can tentatively support three cases:
323 : // - 1 element on the other side, coarse as well (towards less refinement).
324 : // In that case, do not do anything. Two coarse elements sitting next to each other is
325 : // perfect. We can detect this case by looking at the element volumes, with a heuristic
326 : // on the ratio of volumes
327 : // - same number of elements on the other side than the fine elements touching the face
328 : // (refinement was uniform on both sides of the face, we have coarsened one side so far)
329 : // - more elements on the other side than the fine elements touching the face
330 : // (more refinement on that side of the face initially, we are now two levels of
331 : // refinement away)
332 : // TODO: That last case
333 38304 : unsigned int fine_side_index = 0;
334 38304 : const auto coarse_side_center = parent_ptr->side_ptr(side_index)->vertex_average();
335 38304 : Real min_distance = std::numeric_limits<Real>::max();
336 : // There might be a better way to find this index. Smallest distance should work
337 261216 : for (const auto side_index : make_range(fine_el->n_sides()))
338 : {
339 : // only two sides (quad), three sides (hex) also own the coarse node
340 222912 : if (fine_el->side_ptr(side_index)->get_node_index(coarse_node) == libMesh::invalid_uint)
341 111456 : continue;
342 : const auto dist =
343 111456 : (fine_el->side_ptr(side_index)->vertex_average() - coarse_side_center).norm_sq();
344 111456 : if (min_distance > dist)
345 : {
346 63160 : min_distance = dist;
347 63160 : fine_side_index = side_index;
348 : }
349 : }
350 : mooseAssert(min_distance != std::numeric_limits<Real>::max(),
351 : "We should have found a side");
352 :
353 : // We cannot use the neighbor pointer from the fine element, or else wont be able to
354 : // deal with non-conformal meshes that are disjoint at this location
355 : // Instead we offset a little and use a point locator
356 : Point offset_point =
357 76608 : fine_el->side_ptr(fine_side_index)->vertex_average() +
358 0 : 100 * TOLERANCE *
359 114912 : (fine_el->side_ptr(fine_side_index)->vertex_average() - fine_el->vertex_average());
360 38304 : auto pl = mesh_copy->sub_point_locator();
361 38304 : pl->enable_out_of_mesh_mode();
362 38304 : auto const_neighbor = (*pl)(offset_point);
363 38304 : pl->disable_out_of_mesh_mode();
364 :
365 : // We're at a boundary
366 38304 : if (!const_neighbor)
367 4688 : continue;
368 :
369 : // Get a non-const element since it will be a candidate for deletion
370 33616 : auto neighbor_fine_elem = mesh_copy->elem_ptr(const_neighbor->id());
371 :
372 : // Point locator finding a fine element inside
373 33616 : if (fine_elements.find(neighbor_fine_elem) != fine_elements.end())
374 0 : continue;
375 :
376 : // Get the interior node for the next tentative coarse element
377 : // We can just use the index to get it from the next tentative fine element
378 33616 : const auto neighbor_coarse_node_index = neighbor_fine_elem->get_node_index(coarse_node);
379 : // Node is not shared between the coarse element and its fine neighbors.
380 : // The mesh should probably be stitched before attempting coarsening
381 33616 : if (neighbor_coarse_node_index == libMesh::invalid_uint)
382 : {
383 1720 : mooseInfoRepeated("Coarse element node " + Moose::stringify(*coarse_node) +
384 : " does not seem shared with any element other than the coarse element. "
385 : "Is the mesh will stitched? Or are there non-conformalities?");
386 1720 : continue;
387 : }
388 31896 : const auto opposite_node_index = MeshCoarseningUtils::getOppositeNodeIndex(
389 31896 : neighbor_fine_elem->type(), neighbor_coarse_node_index);
390 31896 : auto neighbor_interior_node = neighbor_fine_elem->node_ptr(opposite_node_index);
391 :
392 : // avoid attempting to coarsen again an element we've already coarsened
393 31896 : if (coarse_elems.find(neighbor_fine_elem) == coarse_elems.end())
394 : {
395 : // dont add a candidate if it's too early to coarsen it (will be coarsened on next step)
396 : // dont add a candidate if they are close to the size of a coarsened element already
397 51352 : for (const auto i : index_range(block_ids))
398 49104 : if (block_ids[i] == neighbor_fine_elem->subdomain_id() &&
399 66712 : coarsening[i] > max - coarse_step - 1 &&
400 17608 : std::abs(neighbor_fine_elem->volume()) < std::abs(_max_vol_ratio * fine_el_volume))
401 : {
402 16976 : candidate_pairs.insert(std::make_pair(neighbor_fine_elem, neighbor_interior_node));
403 16976 : break;
404 : }
405 : }
406 38304 : }
407 :
408 : // Delete the elements used to build the coarse element
409 56592 : for (auto & fine_elem : fine_elements)
410 : {
411 49920 : if (!fine_elem)
412 0 : continue;
413 :
414 : // We dont delete nodes in the fine elements as they get removed during renumbering/
415 : // remove_orphaned_nodes calls in preparing for use
416 49920 : mesh_copy->delete_elem(fine_elem);
417 :
418 : // Clean up the list of candidates from any deleted elements
419 870486 : for (auto iter = candidate_pairs.begin(); iter != candidate_pairs.end();)
420 : {
421 820566 : if (iter->first == fine_elem)
422 : {
423 6088 : iter = candidate_pairs.erase(iter);
424 6088 : continue;
425 : }
426 814478 : ++iter;
427 : }
428 : }
429 :
430 : // Contract to remove the elements that were marked for deletion
431 6672 : mesh_copy->contract();
432 : // Prepare for use to refresh the element neighbors
433 6672 : mesh_copy->prepare_for_use();
434 16064 : }
435 :
436 : // We pick the configuration (eg starting node) for which we managed to coarsen the most
437 : // This isn't the best idea, as some coarsening could be invalid (non-conformalities)
438 : // Maybe we should examine for non-conformality here to make a decision?
439 : // It's expensive to do so in a global mesh-wide check though, maybe if we baked that check
440 : // into the coarsening work it would be more reasonable.
441 672 : if (_verbose)
442 576 : _console << "Step " << coarse_step + 1 << " attempt #" << start_node_index << " created "
443 576 : << coarse_elems.size() << " coarse elements." << std::endl;
444 672 : if (int(coarse_elems.size()) > max_num_coarsened)
445 : {
446 184 : mesh_return = std::move(mesh_copy);
447 184 : max_num_coarsened = coarse_elems.size();
448 : }
449 672 : }
450 112 : if (_verbose)
451 96 : _console << "Step " << coarse_step + 1 << " created " << max_num_coarsened
452 96 : << " coarsened elements in its most successful attempt." << std::endl;
453 112 : coarse_step++;
454 112 : return recursiveCoarsen(block_ids, mesh_return, coarsening, max, coarse_step);
455 112 : }
|
