Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/simplex_refiner.h"
21 :
22 : #include "libmesh/boundary_info.h"
23 : #include "libmesh/face_tri3.h"
24 : #include "libmesh/int_range.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/mesh_base.h"
27 : #include "libmesh/mesh_communication.h"
28 : #include "libmesh/mesh_tools.h"
29 : #include "libmesh/node.h"
30 : #include "libmesh/parallel_algebra.h"
31 : #include "libmesh/remote_elem.h"
32 : #include "libmesh/unstructured_mesh.h"
33 : #include "libmesh/utility.h"
34 :
35 : // TIMPI includes
36 : #include "timpi/parallel_implementation.h"
37 : #include "timpi/parallel_sync.h"
38 :
39 :
40 : namespace libMesh
41 : {
42 :
43 : //-----------------------------------------------------------------
44 : // Mesh refinement methods
45 94 : SimplexRefiner::SimplexRefiner (UnstructuredMesh & mesh) :
46 86 : _mesh(mesh),
47 94 : _desired_volume(0.1)
48 94 : {}
49 :
50 :
51 :
52 94 : bool SimplexRefiner::refine_elements()
53 : {
54 94 : if (!_mesh.is_prepared())
55 0 : _mesh.prepare_for_use();
56 :
57 : // We should add more options later: to refine via circumcenter
58 : // insertion instead of edge center insertion, to do Delaunay swaps,
59 : // etc.
60 94 : return this->refine_via_edges();
61 : }
62 :
63 :
64 :
65 44097 : bool SimplexRefiner::should_refine_elem(Elem & elem)
66 : {
67 44097 : const Real min_volume_target = this->desired_volume();
68 44097 : FunctionBase<Real> *volume_func = this->get_desired_volume_function();
69 :
70 : // If this isn't a question, why are we here?
71 3212 : libmesh_assert(min_volume_target > 0 ||
72 : volume_func != nullptr);
73 :
74 : // If we don't have position-dependent volume targets we can make a
75 : // decision quickly
76 44097 : const Real volume = elem.volume();
77 :
78 : if (volume < TOLERANCE*TOLERANCE)
79 : libmesh_warning
80 : ("Warning: trying to refine sliver element with volume " <<
81 : volume << "!\n");
82 :
83 44097 : if (!volume_func)
84 44097 : return (volume > min_volume_target);
85 :
86 : // If we do?
87 : //
88 : // See if we're meeting the local volume target at all the elem
89 : // vertices first
90 0 : for (auto v : make_range(elem.n_vertices()))
91 : {
92 : // If we have an auto volume function, we'll use it and override other volume options
93 0 : const Real local_volume_target = (*volume_func)(elem.point(v));
94 0 : libmesh_error_msg_if
95 : (local_volume_target <= 0,
96 : "Non-positive desired element volumes are unachievable");
97 0 : if (volume > local_volume_target)
98 0 : return true;
99 : }
100 :
101 : // If our vertices are happy, it's still possible that our interior
102 : // isn't. Are we allowed not to bother checking it?
103 0 : if (!min_volume_target)
104 0 : return false;
105 :
106 0 : libmesh_not_implemented_msg
107 : ("Combining a minimum desired_volume with an volume function isn't yet supported.");
108 : }
109 :
110 :
111 44097 : std::size_t SimplexRefiner::refine_via_edges(Elem & elem,
112 : dof_id_type coarse_id)
113 : {
114 : // Do we need to be refined for our own sake, or only if our
115 : // neighbors were refined?
116 44097 : const bool should_refine = this->should_refine_elem(elem);
117 :
118 : // We only currently handle coarse conforming meshes here.
119 : // Uniformly refined meshes should be flattened before using this.
120 3212 : libmesh_assert_equal_to(elem.level(), 0u);
121 :
122 : // Adding Tets or Tri6/Tri7 wouldn't be too much harder. Quads or
123 : // other elements will get tricky.
124 3212 : libmesh_assert_equal_to(elem.type(), TRI3);
125 :
126 44097 : const int n_sides = elem.n_sides();
127 3212 : int side_to_refine = 0;
128 3212 : Real length_to_refine = 0;
129 44097 : std::pair<Node *, Node *> vertices_to_refine;
130 :
131 : // Find the longest side and refine it first. We might change that
132 : // heuristic at some point, to enable anisotropic refinement.
133 176388 : for (int side : make_range(n_sides))
134 : {
135 : std::pair<Node *, Node *> vertices
136 132291 : {elem.node_ptr((side+1)%n_sides),
137 19272 : elem.node_ptr(side)};
138 :
139 132291 : if ((Point&)(*vertices.first) >
140 9636 : (Point&)(*vertices.second))
141 4820 : std::swap(vertices.first, vertices.second);
142 :
143 132291 : const auto sidevec = *(Point*)vertices.first -
144 132291 : *(Point*)vertices.second;
145 122655 : Real length = sidevec.norm();
146 :
147 : // We might need to ask the owner of a coarse ghost element
148 : // about whether it saw any refinement forced by a coarse
149 : // neighbor with a smaller desired area.
150 132291 : if (const processor_id_type pid = elem.processor_id();
151 137109 : pid != _mesh.processor_id() && !_mesh.is_serial() &&
152 137109 : !_desired_volume_func.get() &&
153 74763 : elem.id() == coarse_id)
154 1980 : edge_queries[pid].emplace_back(vertices.first->id(),
155 990 : vertices.second->id());
156 :
157 : // Test should_refine, but also test for any edges that were
158 : // already refined by neighbors, which might happen even if
159 : // !should_refine in cases where we're refining non-uniformly.
160 139719 : if (length > length_to_refine &&
161 99633 : (should_refine ||
162 17064 : new_nodes.find(vertices) != new_nodes.end()))
163 : {
164 1496 : side_to_refine = side;
165 1496 : length_to_refine = length;
166 2992 : vertices_to_refine = vertices;
167 : }
168 : }
169 :
170 : // We might not refine anything.
171 44097 : if (!length_to_refine)
172 2348 : return 0;
173 :
174 : // But if we must, then do it.
175 : //
176 : // Find all the edge neighbors we can. For a triangle there'll just
177 : // be the one, but we want to extend to tets later.
178 1728 : std::vector<Elem *> neighbors;
179 13158 : if (Elem * neigh = elem.neighbor_ptr(side_to_refine); neigh)
180 11310 : neighbors.push_back(neigh);
181 :
182 : Node * midedge_node;
183 12294 : if (auto it = new_nodes.find(vertices_to_refine);
184 864 : it != new_nodes.end())
185 : {
186 5535 : midedge_node = it->second;
187 : }
188 : else
189 : {
190 : // Add with our own processor id first, so we can request a temporary
191 : // id that won't conflict with anyone else's canonical ids later
192 : midedge_node =
193 6759 : _mesh.add_point((*(Point*)vertices_to_refine.first +
194 6759 : *(Point*)vertices_to_refine.second)/2,
195 : DofObject::invalid_id,
196 912 : _mesh.processor_id());
197 :
198 : // Then give the new node the best pid we can.
199 6759 : midedge_node->processor_id() = elem.processor_id();
200 12534 : for (const Elem * neigh : neighbors)
201 5775 : if (neigh != remote_elem)
202 5519 : midedge_node->processor_id() = std::min(midedge_node->processor_id(),
203 5927 : neigh->processor_id());
204 :
205 6759 : new_nodes[vertices_to_refine] = midedge_node;
206 : }
207 :
208 : // Good for Tris and Tets; we'll need more for Quads.
209 864 : constexpr int max_subelems = 2;
210 26316 : std::unique_ptr<Tri3> subelem[max_subelems];
211 :
212 : // We'll switch(elem->type()) here eventually
213 22860 : subelem[0] = std::make_unique<Tri3>();
214 13158 : subelem[0]->set_node(0, elem.node_ptr(side_to_refine));
215 12294 : subelem[0]->set_node(1, midedge_node);
216 13158 : subelem[0]->set_node(2, elem.node_ptr((side_to_refine+2)%3));
217 :
218 22860 : subelem[1] = std::make_unique<Tri3>();
219 13158 : subelem[1]->set_node(0, elem.node_ptr((side_to_refine+2)%3));
220 12294 : subelem[1]->set_node(1, midedge_node);
221 13158 : subelem[1]->set_node(2, elem.node_ptr((side_to_refine+1)%3));
222 :
223 : // Preserve boundary and sort-of-preserve neighbor information. We
224 : // need remote_elem neighbors to avoid breaking distributed meshes
225 : // (which can't reconstruct them from scratch), and we need at least
226 : // placeholder neighbors to make sure we'll have good processor_id()
227 : // options for new nodes in future splits of the same edge.
228 1728 : subelem[0]->set_neighbor(1, proxy_elements[elem.processor_id()].get());
229 1728 : subelem[1]->set_neighbor(0, proxy_elements[elem.processor_id()].get());
230 :
231 12294 : BoundaryInfo & boundary_info = _mesh.get_boundary_info();
232 864 : std::vector<boundary_id_type> bc_ids;
233 12294 : boundary_info.boundary_ids(&elem, side_to_refine, bc_ids);
234 12294 : boundary_info.add_side(subelem[0].get(), 0, bc_ids);
235 12294 : boundary_info.add_side(subelem[1].get(), 1, bc_ids);
236 2592 : subelem[0]->set_neighbor(0, elem.neighbor_ptr(side_to_refine));
237 1728 : subelem[1]->set_neighbor(1, elem.neighbor_ptr(side_to_refine));
238 :
239 12294 : boundary_info.boundary_ids(&elem, (side_to_refine+1)%3, bc_ids);
240 12294 : boundary_info.add_side(subelem[1].get(), 2, bc_ids);
241 2592 : subelem[1]->set_neighbor(2, elem.neighbor_ptr((side_to_refine+1)%3));
242 :
243 12294 : boundary_info.boundary_ids(&elem, (side_to_refine+2)%3, bc_ids);
244 12294 : boundary_info.add_side(subelem[0].get(), 2, bc_ids);
245 2592 : subelem[0]->set_neighbor(2, elem.neighbor_ptr((side_to_refine+2)%3));
246 :
247 : // Be sure the correct data is set for all subelems.
248 12294 : const unsigned int nei = elem.n_extra_integers();
249 36882 : for (unsigned int i=0; i != max_subelems; ++i)
250 26316 : if (subelem[i])
251 : {
252 : // We'd love to assert that we don't have a sliver here,
253 : // but maybe our input had a sliver and we can't fix it.
254 1728 : libmesh_assert_greater_equal(subelem[i]->volume() + TOLERANCE,
255 : elem.volume() / 2);
256 :
257 : // Copy any extra element data. Since the subelements
258 : // haven't been added to the mesh yet any allocation has
259 : // to be done manually.
260 24588 : subelem[i]->add_extra_integers(nei);
261 24588 : for (unsigned int ei=0; ei != nei; ++ei)
262 0 : subelem[ei]->set_extra_integer(ei, elem.get_extra_integer(ei));
263 :
264 24588 : subelem[i]->inherit_data_from(elem);
265 : }
266 :
267 12294 : boundary_info.remove(&elem);
268 :
269 12294 : if (auto it = this->added_elements.find(coarse_id);
270 864 : it != this->added_elements.end())
271 : {
272 832 : auto & added = it->second;
273 22680 : if (auto sub_it = added.find(elem.vertex_average());
274 832 : sub_it != added.end())
275 : {
276 11756 : if (&elem == sub_it->second)
277 832 : added.erase(sub_it);
278 : }
279 : }
280 :
281 : // We only add an element to new_elements *right* before
282 : // checking it for further refinement, so we only need to check the
283 : // back() to see where the elem here came from.
284 13152 : if (!this->new_elements.empty() &&
285 858 : &elem == this->new_elements.back().get())
286 : {
287 748 : this->new_elements.pop_back();
288 10855 : this->coarse_parent.erase(&elem);
289 : }
290 : else
291 1439 : _mesh.delete_elem(&elem);
292 :
293 : // We refined at least ourselves
294 864 : std::size_t n_refined_elements = 1;
295 :
296 36882 : for (unsigned int i=0; i != max_subelems; ++i)
297 : {
298 26316 : Elem * add_elem = subelem[i].get();
299 24588 : added_elements[coarse_id].emplace(add_elem->vertex_average(),
300 1728 : add_elem);
301 24588 : this->new_elements.push_back(std::move(subelem[i]));
302 24588 : this->coarse_parent[add_elem] = coarse_id;
303 24588 : n_refined_elements +=
304 24588 : this->refine_via_edges(*add_elem, coarse_id);
305 : }
306 :
307 864 : return n_refined_elements;
308 31698 : }
309 :
310 :
311 :
312 94 : std::size_t SimplexRefiner::refine_via_edges()
313 : {
314 4 : this->new_nodes.clear();
315 90 : this->new_elements.clear();
316 :
317 : // We need to look at neighbors' pids to determine new node pids,
318 : // but we might be deleting elements and leaving dangling neighbor
319 : // pointers. Put proxy neighbors in those pointers places instead.
320 1196 : for (auto & elem : _mesh.element_ptr_range())
321 2152 : for (auto n : make_range(elem->n_neighbors()))
322 : {
323 1614 : Elem * neigh = elem->neighbor_ptr(n);
324 1614 : if (!neigh || neigh == remote_elem)
325 408 : continue;
326 :
327 1206 : const processor_id_type p = neigh->processor_id();
328 80 : std::unique_ptr<Elem> & proxy_neighbor = proxy_elements[p];
329 1206 : if (!proxy_neighbor.get())
330 748 : proxy_neighbor = Elem::build(NODEELEM);
331 :
332 1286 : elem->set_neighbor(n, proxy_neighbor.get());
333 86 : }
334 :
335 : // If we have a _desired_volume_func then we might still have
336 : // some unsplit edges that could have been split by their
337 : // remote_elem neighbors, and we'll need to query those.
338 : auto edge_gather_functor =
339 352 : [this]
340 : (processor_id_type,
341 : const std::vector<std::pair<dof_id_type, dof_id_type>> & edges,
342 1254 : std::vector<refinement_datum> & edge_refinements)
343 : {
344 : // Fill those requests
345 0 : const std::size_t query_size = edges.size();
346 352 : edge_refinements.resize(query_size);
347 1606 : for (std::size_t i=0; i != query_size; ++i)
348 : {
349 1254 : auto vertex_ids = edges[i];
350 : std::pair<Node *, Node *> vertices
351 1254 : {_mesh.node_ptr(vertex_ids.first),
352 1254 : _mesh.node_ptr(vertex_ids.second)};
353 1254 : fill_refinement_datum(vertices, edge_refinements[i]);
354 : }
355 442 : };
356 :
357 : auto edge_action_functor =
358 352 : [this]
359 : (processor_id_type,
360 : const std::vector<std::pair<dof_id_type, dof_id_type>> &,
361 : const std::vector<refinement_datum> & edge_refinements
362 4466 : )
363 : {
364 1606 : for (const auto & one_edges_refinements : edge_refinements)
365 5720 : for (const auto & refinement : one_edges_refinements)
366 : {
367 : std::pair<Node *, Node *> vertices
368 4466 : {_mesh.node_ptr(std::get<0>(refinement)),
369 4466 : _mesh.node_ptr(std::get<1>(refinement))};
370 :
371 4466 : if ((Point&)(*vertices.first) >
372 0 : (Point&)(*vertices.second))
373 0 : std::swap(vertices.first, vertices.second);
374 :
375 4466 : if (auto it = new_nodes.find(vertices);
376 0 : it != new_nodes.end())
377 : {
378 4450 : _mesh.renumber_node(it->second->id(),
379 0 : std::get<2>(refinement));
380 4450 : it->second->processor_id() = std::get<3>(refinement);
381 : }
382 : else
383 : {
384 : Node * new_node =
385 32 : _mesh.add_point((*(Point*)vertices.first +
386 16 : *(Point*)vertices.second)/2,
387 0 : std::get<2>(refinement),
388 0 : std::get<3>(refinement));
389 16 : new_nodes[vertices] = new_node;
390 : }
391 : }
392 90 : };
393 :
394 4 : std::size_t refined_elements = 0;
395 4 : std::size_t newly_refined_elements = 0;
396 6 : do
397 : {
398 211 : newly_refined_elements = 0;
399 :
400 : // Refinement results should agree on ghost elements, except for
401 : // id and unique_id assignment
402 19921 : for (auto & elem : _mesh.element_ptr_range())
403 : {
404 1484 : auto it = coarse_parent.find(elem);
405 : const dof_id_type coarse_id =
406 20961 : (it == coarse_parent.end()) ?
407 1990 : elem->id() : it->second;
408 :
409 19509 : newly_refined_elements +=
410 19509 : this->refine_via_edges(*elem, coarse_id);
411 191 : }
412 211 : refined_elements += newly_refined_elements;
413 13944 : for (auto & elem : this->new_elements)
414 15693 : _mesh.add_elem(std::move(elem));
415 10 : this->new_elements.clear();
416 :
417 : // Yeah, this is just a lower bound in parallel
418 211 : _mesh.comm().max(newly_refined_elements);
419 :
420 211 : if (newly_refined_elements && !_mesh.is_replicated())
421 : {
422 96 : bool have_edge_queries = !edge_queries.empty();
423 96 : _mesh.comm().max(have_edge_queries);
424 0 : refinement_datum * refinement_data_ex = nullptr;
425 96 : if (have_edge_queries)
426 : Parallel::pull_parallel_vector_data
427 93 : (_mesh.comm(), edge_queries, edge_gather_functor,
428 : edge_action_functor, refinement_data_ex);
429 : }
430 : }
431 211 : while (newly_refined_elements);
432 :
433 : // If we're replicated then we should be done here.
434 94 : if (!_mesh.is_replicated())
435 : {
436 : // If we're not replicated then we might have a bunch of nodes
437 : // and elements that still have temporary id and unique_id
438 : // values. Give them permanent ones.
439 0 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> elems_to_query;
440 506 : for (const auto & [coarse_id, added_elem_map] : added_elements)
441 : {
442 426 : if (added_elem_map.empty())
443 96 : continue;
444 :
445 : const processor_id_type pid =
446 426 : added_elem_map.begin()->second->processor_id();
447 426 : if (pid == _mesh.processor_id())
448 96 : continue;
449 :
450 330 : elems_to_query[pid].push_back(coarse_id);
451 : }
452 :
453 : // Return fine element data based on vertex_average()
454 : typedef
455 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
456 : std::vector<std::tuple<Point, dof_id_type, unique_id_type>>
457 : #else
458 : std::vector<std::tuple<Point, dof_id_type>>
459 : #endif
460 : elem_refinement_datum;
461 :
462 : auto added_gather_functor =
463 288 : [this]
464 : (processor_id_type,
465 : const std::vector<dof_id_type> & coarse_elems,
466 330 : std::vector<elem_refinement_datum> & coarse_refinements)
467 : {
468 : // Fill those requests
469 0 : const std::size_t query_size = coarse_elems.size();
470 288 : coarse_refinements.resize(query_size);
471 618 : for (auto i : make_range(query_size))
472 : {
473 330 : const dof_id_type coarse_id = coarse_elems[i];
474 : const auto & added =
475 330 : libmesh_map_find(added_elements, coarse_id);
476 :
477 7722 : for (auto [vertex_avg, elem] : added)
478 : {
479 0 : coarse_refinements[i].emplace_back
480 7392 : (vertex_avg,
481 7392 : elem->id()
482 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
483 7392 : , elem->unique_id()
484 : #endif
485 0 : );
486 : }
487 : }
488 288 : };
489 :
490 : auto added_action_functor =
491 288 : [this]
492 : (processor_id_type,
493 : const std::vector<dof_id_type> & coarse_elems,
494 7722 : const std::vector<elem_refinement_datum> & coarse_refinements)
495 : {
496 0 : const std::size_t query_size = coarse_elems.size();
497 618 : for (auto i : make_range(query_size))
498 : {
499 330 : const dof_id_type coarse_id = coarse_elems[i];
500 0 : const auto & refinement_data = coarse_refinements[i];
501 : const auto & our_added =
502 330 : libmesh_map_find(added_elements, coarse_id);
503 330 : for (auto [vertex_avg, id
504 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
505 : , unique_id
506 : #endif
507 7722 : ] : refinement_data)
508 : {
509 7392 : Elem & our_elem = *libmesh_map_find(our_added,
510 : vertex_avg);
511 0 : libmesh_assert_equal_to(our_elem.vertex_average(),
512 : vertex_avg);
513 :
514 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
515 7392 : our_elem.set_unique_id(unique_id);
516 : #endif
517 7392 : _mesh.renumber_elem(our_elem.id(), id);
518 : }
519 : }
520 288 : };
521 :
522 0 : elem_refinement_datum * refinement_data_ex = nullptr;
523 : Parallel::pull_parallel_vector_data
524 80 : (_mesh.comm(), elems_to_query, added_gather_functor,
525 : added_action_functor, refinement_data_ex);
526 :
527 : // That took care of our element ids; now get node ids.
528 : MeshCommunication mc;
529 80 : mc.make_node_proc_ids_parallel_consistent (_mesh);
530 80 : mc.make_node_ids_parallel_consistent (_mesh);
531 80 : mc.make_node_unique_ids_parallel_consistent (_mesh);
532 : }
533 :
534 : // In theory the operations on ghost elements are, after remote edge
535 : // splittings are known, embarrassingly parallel. In practice ...
536 : // let's verify that we haven't just desynchronized our mesh.
537 : #ifdef DEBUG
538 4 : MeshTools::libmesh_assert_equal_points(_mesh);
539 4 : MeshTools::libmesh_assert_equal_connectivity (_mesh);
540 : # ifdef LIBMESH_ENABLE_UNIQUE_ID
541 4 : MeshTools::libmesh_assert_valid_unique_ids(_mesh);
542 : # endif
543 : #endif
544 :
545 94 : return refined_elements;
546 : }
547 :
548 :
549 :
550 0 : void SimplexRefiner::set_desired_volume_function
551 : (FunctionBase<Real> * desired)
552 : {
553 0 : if (desired)
554 0 : _desired_volume_func = desired->clone();
555 : else
556 0 : _desired_volume_func.reset();
557 0 : }
558 :
559 :
560 :
561 44097 : FunctionBase<Real> * SimplexRefiner::get_desired_volume_function ()
562 : {
563 44097 : return _desired_volume_func.get();
564 : }
565 :
566 :
567 10186 : void SimplexRefiner::fill_refinement_datum(std::pair<Node *, Node *> vertices,
568 : refinement_datum & vec)
569 : {
570 0 : auto it = new_nodes.find(vertices);
571 10186 : if (it == new_nodes.end())
572 0 : return;
573 :
574 8932 : vec.emplace_back(vertices.first->id(),
575 8932 : vertices.second->id(),
576 4466 : it->second->id(),
577 4466 : it->second->processor_id());
578 :
579 0 : std::pair<Node *, Node *> subedge {vertices.first, it->second};
580 4466 : if ((Point&)(*subedge.first) >
581 0 : (Point&)(*subedge.second))
582 0 : std::swap(subedge.first, subedge.second);
583 4466 : fill_refinement_datum(subedge, vec);
584 :
585 0 : subedge = std::make_pair(it->second, vertices.second);
586 4466 : if ((Point&)(*subedge.first) >
587 0 : (Point&)(*subedge.second))
588 0 : std::swap(subedge.first, subedge.second);
589 4466 : fill_refinement_datum(subedge, vec);
590 : }
591 :
592 :
593 : } // namespace libMesh
594 :
|
