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 :
20 : // C++ includes
21 :
22 : // Local includes
23 : #include "libmesh/mesh_tools.h"
24 : #include "libmesh/mesh_subdivision_support.h"
25 : #include "libmesh/boundary_info.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 37177 : void MeshTools::Subdivision::find_one_ring(const Tri3Subdivision * elem,
32 : std::vector<const Node *> & nodes)
33 : {
34 3104 : libmesh_assert(elem->is_subdivision_updated());
35 3104 : libmesh_assert(elem->get_ordered_node(0));
36 :
37 37177 : unsigned int valence = elem->get_ordered_valence(0);
38 37177 : nodes.resize(valence + 6);
39 :
40 : // The first three vertices in the patch are the ones from the element triangle
41 37177 : nodes[0] = elem->get_ordered_node(0);
42 37177 : nodes[1] = elem->get_ordered_node(1);
43 37177 : nodes[valence] = elem->get_ordered_node(2);
44 :
45 37177 : const unsigned int nn0 = elem->local_node_number(nodes[0]->id());
46 :
47 40281 : const Tri3Subdivision * nb = dynamic_cast<const Tri3Subdivision *>(elem->neighbor_ptr(nn0));
48 3104 : libmesh_assert(nb);
49 :
50 3104 : unsigned int j, i = 1;
51 :
52 12312 : do
53 : {
54 184715 : ++i;
55 184715 : j = nb->local_node_number(nodes[0]->id());
56 200131 : nodes[i] = nb->node_ptr(next[j]);
57 30832 : nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
58 184715 : } while (nb != elem);
59 :
60 : /* for nodes connected with N (= valence[0]) */
61 37177 : nb = static_cast<const Tri3Subdivision *>(elem->neighbor_ptr(next[nn0]));
62 37177 : j = nb->local_node_number(nodes[1]->id());
63 40281 : nodes[valence+1] = nb->node_ptr(next[j]);
64 :
65 6208 : nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(next[j]));
66 37177 : j = nb->local_node_number(nodes[valence+1]->id());
67 40281 : nodes[valence+4] = nb->node_ptr(next[j]);
68 :
69 6208 : nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(next[j]));
70 37177 : j = nb->local_node_number(nodes[valence+4]->id());
71 40281 : nodes[valence+5] = nb->node_ptr(next[j]);
72 :
73 : /* for nodes connected with 1 */
74 6208 : nb = static_cast<const Tri3Subdivision *>(elem->neighbor_ptr(next[nn0]));
75 37177 : j = nb->local_node_number(nodes[1]->id());
76 : // nodes[valence+1] has been determined already
77 :
78 6208 : nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
79 37177 : j = nb->local_node_number(nodes[1]->id());
80 40281 : nodes[valence+2] = nb->node_ptr(next[j]);
81 :
82 6208 : nb = static_cast<const Tri3Subdivision *>(nb->neighbor_ptr(j));
83 37177 : j = nb->local_node_number(nodes[1]->id());
84 43385 : nodes[valence+3] = nb->node_ptr(next[j]);
85 :
86 40281 : return;
87 : }
88 :
89 :
90 71 : void MeshTools::Subdivision::all_subdivision(MeshBase & mesh)
91 : {
92 : const bool mesh_has_boundary_data =
93 2 : (mesh.get_boundary_info().n_boundary_ids() > 0);
94 :
95 4 : std::vector<Elem *> new_boundary_elements;
96 4 : std::vector<short int> new_boundary_sides;
97 4 : std::vector<boundary_id_type> new_boundary_ids;
98 :
99 : // Container to catch ids handed back from BoundaryInfo
100 4 : std::vector<boundary_id_type> ids;
101 :
102 35468 : for (const auto & elem : mesh.element_ptr_range())
103 : {
104 512 : libmesh_assert_equal_to(elem->type(), TRI3);
105 :
106 18176 : auto tri = Elem::build_with_id(TRI3SUBDIVISION, elem->id());
107 18176 : tri->subdomain_id() = elem->subdomain_id();
108 18688 : tri->set_node(0, elem->node_ptr(0));
109 18688 : tri->set_node(1, elem->node_ptr(1));
110 18688 : tri->set_node(2, elem->node_ptr(2));
111 :
112 18176 : if (mesh_has_boundary_data)
113 : {
114 0 : for (auto side : elem->side_index_range())
115 : {
116 0 : mesh.get_boundary_info().boundary_ids(elem, side, ids);
117 :
118 0 : for (const auto & id : ids)
119 : {
120 : // add the boundary id to the list of new boundary ids
121 0 : new_boundary_ids.push_back(id);
122 0 : new_boundary_elements.push_back(tri.get());
123 0 : new_boundary_sides.push_back(side);
124 : }
125 : }
126 :
127 : // remove the original element from the BoundaryInfo structure
128 0 : mesh.get_boundary_info().remove(elem);
129 : }
130 :
131 19200 : mesh.insert_elem(std::move(tri));
132 17219 : }
133 71 : mesh.prepare_for_use();
134 :
135 71 : if (mesh_has_boundary_data)
136 : {
137 : // If the old mesh had boundary data, the new mesh better have some too.
138 0 : libmesh_assert_greater(new_boundary_elements.size(), 0);
139 :
140 : // We should also be sure that the lengths of the new boundary data vectors
141 : // are all the same.
142 0 : libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_elements.size());
143 0 : libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_ids.size());
144 :
145 : // Add the new boundary info to the mesh.
146 0 : for (auto s : index_range(new_boundary_elements))
147 0 : mesh.get_boundary_info().add_side(new_boundary_elements[s],
148 0 : new_boundary_sides[s],
149 0 : new_boundary_ids[s]);
150 : }
151 :
152 71 : mesh.prepare_for_use();
153 71 : }
154 :
155 :
156 71 : void MeshTools::Subdivision::prepare_subdivision_mesh(MeshBase & mesh, bool ghosted)
157 : {
158 71 : mesh.prepare_for_use();
159 :
160 : // convert all mesh elements to subdivision elements
161 71 : all_subdivision(mesh);
162 :
163 71 : if (!ghosted)
164 : {
165 : // add the ghost elements for the boundaries
166 71 : add_boundary_ghosts(mesh);
167 : }
168 : else
169 : {
170 : // This assumes that the mesh already has the ghosts. Only tagging them is required here.
171 0 : tag_boundary_ghosts(mesh);
172 : }
173 :
174 71 : mesh.prepare_for_use();
175 :
176 4 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
177 71 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
178 :
179 : // compute the node valences
180 25118 : for (auto & node : mesh.node_ptr_range())
181 : {
182 724 : std::vector<const Node *> neighbors;
183 12851 : MeshTools::find_nodal_neighbors(mesh, *node, nodes_to_elem_map, neighbors);
184 : const unsigned int valence =
185 724 : cast_int<unsigned int>(neighbors.size());
186 362 : libmesh_assert_greater(valence, 1);
187 12851 : node->set_valence(valence);
188 67 : }
189 :
190 44852 : for (auto & elem : mesh.element_ptr_range())
191 : {
192 23004 : Tri3Subdivision * tri3s = dynamic_cast<Tri3Subdivision *>(elem);
193 648 : libmesh_assert(tri3s);
194 23004 : if (!tri3s->is_ghost())
195 18176 : tri3s->prepare_subdivision_properties();
196 67 : }
197 71 : }
198 :
199 :
200 0 : void MeshTools::Subdivision::tag_boundary_ghosts(MeshBase & mesh)
201 : {
202 0 : for (auto & elem : mesh.element_ptr_range())
203 : {
204 0 : libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
205 :
206 0 : Tri3Subdivision * sd_elem = static_cast<Tri3Subdivision *>(elem);
207 0 : for (auto i : elem->side_index_range())
208 : {
209 0 : if (elem->neighbor_ptr(i) == nullptr)
210 : {
211 0 : sd_elem->set_ghost(true);
212 : // set all other neighbors to ghosts as well
213 0 : if (elem->neighbor_ptr(next[i]))
214 : {
215 0 : Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(next[i]));
216 0 : nb->set_ghost(true);
217 : }
218 0 : if (elem->neighbor_ptr(prev[i]))
219 : {
220 0 : Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(prev[i]));
221 0 : nb->set_ghost(true);
222 : }
223 : }
224 : }
225 0 : }
226 0 : }
227 :
228 :
229 71 : void MeshTools::Subdivision::add_boundary_ghosts(MeshBase & mesh)
230 : {
231 : static const Real tol = 1e-5;
232 :
233 : // add the mirrored ghost elements (without using iterators, because the mesh is modified in the course)
234 4 : std::vector<Tri3Subdivision *> ghost_elems;
235 4 : std::vector<Node *> ghost_nodes;
236 71 : const unsigned int n_elem = mesh.n_elem();
237 18247 : for (unsigned int eid = 0; eid < n_elem; ++eid)
238 : {
239 18176 : Elem * elem = mesh.elem_ptr(eid);
240 512 : libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
241 :
242 : // If the triangle happens to be in a corner (two boundary
243 : // edges), we perform a counter-clockwise loop by mirroring the
244 : // previous triangle until we come back to the original
245 : // triangle. This prevents degenerated triangles in the mesh
246 : // corners and guarantees that the node in the middle of the
247 : // loop is of valence=6.
248 72704 : for (auto i : elem->side_index_range())
249 : {
250 1536 : libmesh_assert_not_equal_to(elem->neighbor_ptr(i), elem);
251 :
252 56128 : if (elem->neighbor_ptr(i) == nullptr &&
253 2272 : elem->neighbor_ptr(next[i]) == nullptr)
254 : {
255 0 : Elem * nelem = elem;
256 0 : unsigned int k = i;
257 0 : for (unsigned int l=0;l<4;l++)
258 : {
259 : // this is the vertex to be mirrored
260 0 : Point point = nelem->point(k) + nelem->point(next[k]) - nelem->point(prev[k]);
261 :
262 : // Check if the proposed vertex doesn't coincide
263 : // with one of the existing vertices. This is
264 : // necessary because for some triangulations, it can
265 : // happen that two mirrored ghost vertices coincide,
266 : // which would then lead to a zero size ghost
267 : // element below.
268 0 : Node * node = nullptr;
269 0 : for (auto & ghost_node : ghost_nodes)
270 0 : if ((*ghost_node - point).norm() < tol * (elem->point(k) - point).norm())
271 : {
272 0 : node = ghost_node;
273 0 : break;
274 : }
275 :
276 : // add the new vertex only if no other is nearby
277 0 : if (node == nullptr)
278 : {
279 0 : node = mesh.add_point(point);
280 0 : ghost_nodes.push_back(node);
281 : }
282 :
283 0 : auto uelem = Elem::build(TRI3SUBDIVISION);
284 0 : auto newelem = cast_ptr<Tri3Subdivision *>(uelem.get());
285 :
286 : // add the first new ghost element to the list just as in the non-corner case
287 0 : if (l == 0)
288 0 : ghost_elems.push_back(newelem);
289 :
290 0 : newelem->set_node(0, nelem->node_ptr(next[k]));
291 0 : newelem->set_node(1, nelem->node_ptr(k));
292 0 : newelem->set_node(2, node);
293 0 : newelem->set_neighbor(0, nelem);
294 0 : newelem->set_ghost(true);
295 0 : if (l>0)
296 0 : newelem->set_neighbor(2, nullptr);
297 0 : nelem->set_neighbor(k, newelem);
298 :
299 0 : Elem * added_elem = mesh.add_elem(std::move(uelem));
300 0 : mesh.get_boundary_info().add_node(nelem->node_ptr(k), 1);
301 0 : mesh.get_boundary_info().add_node(nelem->node_ptr(next[k]), 1);
302 0 : mesh.get_boundary_info().add_node(nelem->node_ptr(prev[k]), 1);
303 0 : mesh.get_boundary_info().add_node(node, 1);
304 :
305 0 : nelem = added_elem;
306 0 : k = 2 ;
307 0 : }
308 :
309 0 : auto uelem = Elem::build(TRI3SUBDIVISION);
310 0 : auto newelem = cast_ptr<Tri3Subdivision *>(uelem.get());
311 :
312 0 : newelem->set_node(0, elem->node_ptr(next[i]));
313 0 : newelem->set_node(1, nelem->node_ptr(2));
314 0 : newelem->set_node(2, elem->node_ptr(prev[i]));
315 0 : newelem->set_neighbor(0, nelem);
316 0 : nelem->set_neighbor(2, newelem);
317 0 : newelem->set_ghost(true);
318 0 : newelem->set_neighbor(2, elem);
319 0 : elem->set_neighbor(next[i],newelem);
320 :
321 0 : mesh.add_elem(std::move(uelem));
322 :
323 0 : break;
324 0 : }
325 : }
326 :
327 72704 : for (auto i : elem->side_index_range())
328 : {
329 1536 : libmesh_assert_not_equal_to(elem->neighbor_ptr(i), elem);
330 56064 : if (elem->neighbor_ptr(i) == nullptr)
331 : {
332 : // this is the vertex to be mirrored
333 2336 : Point point = elem->point(i) + elem->point(next[i]) - elem->point(prev[i]);
334 :
335 : // Check if the proposed vertex doesn't coincide with
336 : // one of the existing vertices. This is necessary
337 : // because for some triangulations, it can happen that
338 : // two mirrored ghost vertices coincide, which would
339 : // then lead to a zero size ghost element below.
340 2272 : Node * node = nullptr;
341 37488 : for (auto & ghost_node : ghost_nodes)
342 37200 : if ((*ghost_node - point).norm() < tol * (elem->point(i) - point).norm())
343 : {
344 0 : node = ghost_node;
345 0 : break;
346 : }
347 :
348 : // add the new vertex only if no other is nearby
349 2272 : if (node == nullptr)
350 : {
351 2272 : node = mesh.add_point(point);
352 2272 : ghost_nodes.push_back(node);
353 : }
354 :
355 2336 : auto uelem = Elem::build(TRI3SUBDIVISION);
356 2272 : auto newelem = cast_ptr<Tri3Subdivision *>(uelem.get());
357 :
358 2272 : ghost_elems.push_back(newelem);
359 :
360 192 : newelem->set_node(0, elem->node_ptr(next[i]));
361 192 : newelem->set_node(1, elem->node_ptr(i));
362 2272 : newelem->set_node(2, node);
363 128 : newelem->set_neighbor(0, elem);
364 64 : newelem->set_ghost(true);
365 128 : elem->set_neighbor(i, newelem);
366 :
367 2400 : mesh.add_elem(std::move(uelem));
368 2336 : mesh.get_boundary_info().add_node(elem->node_ptr(i), 1);
369 2336 : mesh.get_boundary_info().add_node(elem->node_ptr(next[i]), 1);
370 2336 : mesh.get_boundary_info().add_node(elem->node_ptr(prev[i]), 1);
371 2272 : mesh.get_boundary_info().add_node(node, 1);
372 2144 : }
373 : }
374 : }
375 :
376 : // add the missing ghost elements (connecting new ghost nodes)
377 6 : std::vector<std::unique_ptr<Elem>> missing_ghost_elems;
378 2343 : for (auto & elem : ghost_elems)
379 : {
380 64 : libmesh_assert(elem->is_ghost());
381 :
382 4608 : for (auto i : elem->side_index_range())
383 : {
384 4736 : if (elem->neighbor_ptr(i) == nullptr &&
385 2272 : elem->neighbor_ptr(prev[i]) != nullptr)
386 : {
387 : // go around counter-clockwise
388 64 : Tri3Subdivision * nb1 = static_cast<Tri3Subdivision *>(elem->neighbor_ptr(prev[i]));
389 64 : Tri3Subdivision * nb2 = nb1;
390 64 : unsigned int j = i;
391 64 : unsigned int n_nb = 0;
392 11076 : while (nb1 != nullptr && nb1->id() != elem->id())
393 : {
394 9052 : j = nb1->local_node_number(elem->node_id(i));
395 248 : nb2 = nb1;
396 8804 : nb1 = static_cast<Tri3Subdivision *>(nb1->neighbor_ptr(prev[j]));
397 248 : libmesh_assert(nb1 == nullptr || nb1->id() != nb2->id());
398 8804 : n_nb++;
399 : }
400 :
401 64 : libmesh_assert_not_equal_to(nb2->id(), elem->id());
402 :
403 : // Above, we merged coinciding ghost vertices. Therefore, we need
404 : // to exclude the case where there is no ghost element to add between
405 : // these two (identical) ghost nodes.
406 2400 : if (elem->node_ptr(next[i])->id() == nb2->node_ptr(prev[j])->id())
407 0 : break;
408 :
409 : // If the number of already present neighbors is less than 4, we add another extra element
410 : // so that the node in the middle of the loop ends up being of valence=6.
411 : // This case usually happens when the middle node corresponds to a corner of the original mesh,
412 : // and the extra element below prevents degenerated triangles in the mesh corners.
413 2272 : if (n_nb < 4)
414 : {
415 : // this is the vertex to be mirrored
416 284 : Point point = nb2->point(j) + nb2->point(prev[j]) - nb2->point(next[j]);
417 :
418 : // Check if the proposed vertex doesn't coincide with one of the existing vertices.
419 : // This is necessary because for some triangulations, it can happen that two mirrored
420 : // ghost vertices coincide, which would then lead to a zero size ghost element below.
421 284 : Node * node = nullptr;
422 9798 : for (auto & ghost_node : ghost_nodes)
423 10050 : if ((*ghost_node - point).norm() < tol * (nb2->point(j) - point).norm())
424 : {
425 0 : node = ghost_node;
426 0 : break;
427 : }
428 :
429 : // add the new vertex only if no other is nearby
430 284 : if (node == nullptr)
431 : {
432 284 : node = mesh.add_point(point);
433 284 : ghost_nodes.push_back(node);
434 : }
435 :
436 284 : auto uelem = Elem::build(TRI3SUBDIVISION);
437 8 : auto newelem = cast_ptr<Tri3Subdivision *>(uelem.get());
438 :
439 24 : newelem->set_node(0, nb2->node_ptr(j));
440 24 : newelem->set_node(1, nb2->node_ptr(prev[j]));
441 284 : newelem->set_node(2, node);
442 16 : newelem->set_neighbor(0, nb2);
443 8 : newelem->set_neighbor(1, nullptr);
444 8 : newelem->set_ghost(true);
445 16 : nb2->set_neighbor(prev[j], newelem);
446 :
447 300 : Elem * added_elem = mesh.add_elem(std::move(uelem));
448 292 : mesh.get_boundary_info().add_node(nb2->node_ptr(j), 1);
449 292 : mesh.get_boundary_info().add_node(nb2->node_ptr(prev[j]), 1);
450 284 : mesh.get_boundary_info().add_node(node, 1);
451 :
452 8 : nb2 = cast_ptr<Tri3Subdivision *>(added_elem);
453 292 : j = nb2->local_node_number(elem->node_id(i));
454 268 : }
455 :
456 2336 : auto uelem = Elem::build(TRI3SUBDIVISION);
457 64 : auto newelem = cast_ptr<Tri3Subdivision *>(uelem.get());
458 :
459 2336 : newelem->set_node(0, elem->node_ptr(next[i]));
460 192 : newelem->set_node(1, elem->node_ptr(i));
461 2336 : newelem->set_node(2, nb2->node_ptr(prev[j]));
462 128 : newelem->set_neighbor(0, elem);
463 64 : newelem->set_neighbor(1, nb2);
464 64 : newelem->set_neighbor(2, nullptr);
465 64 : newelem->set_ghost(true);
466 :
467 128 : elem->set_neighbor(i, newelem);
468 128 : nb2->set_neighbor(prev[j], newelem);
469 :
470 64 : missing_ghost_elems.push_back(std::move(uelem));
471 64 : break;
472 2144 : }
473 : } // end side loop
474 : } // end ghost element loop
475 :
476 : // add the missing ghost elements to the mesh
477 2343 : for (auto & elem : missing_ghost_elems)
478 2400 : mesh.add_elem(std::move(elem));
479 138 : }
480 :
481 : } // namespace libMesh
|
