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