Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : #include "libmesh/libmesh_config.h"
19 :
20 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
21 :
22 : // C++ includes
23 :
24 : // Local includes
25 : #include "libmesh/inf_elem_builder.h"
26 : #include "libmesh/libmesh_logging.h"
27 : #include "libmesh/mesh_tools.h"
28 : #include "libmesh/face_inf_quad4.h"
29 : #include "libmesh/face_inf_quad6.h"
30 : #include "libmesh/cell_inf_prism6.h"
31 : #include "libmesh/cell_inf_prism12.h"
32 : #include "libmesh/cell_inf_hex8.h"
33 : #include "libmesh/cell_inf_hex16.h"
34 : #include "libmesh/cell_inf_hex18.h"
35 : #include "libmesh/mesh_base.h"
36 : #include "libmesh/remote_elem.h"
37 :
38 : #include "timpi/parallel_implementation.h"
39 :
40 : namespace libMesh
41 : {
42 :
43 76 : const Point InfElemBuilder::build_inf_elem(bool be_verbose)
44 : {
45 : // determine origin automatically,
46 : // works only if the mesh has no symmetry planes.
47 76 : const BoundingBox b_box = MeshTools::create_bounding_box(_mesh);
48 30 : Point origin = (b_box.first + b_box.second) / 2;
49 :
50 30 : if (be_verbose && _mesh.processor_id() == 0)
51 : {
52 : #ifdef DEBUG
53 2 : libMesh::out << " Determined origin for Infinite Elements:"
54 2 : << std::endl
55 2 : << " ";
56 2 : origin.write_unformatted(libMesh::out);
57 2 : libMesh::out << std::endl;
58 : #endif
59 : }
60 :
61 : // Call the protected implementation function with the
62 : // automatically determined origin.
63 76 : this->build_inf_elem(origin, false, false, false, be_verbose);
64 :
65 : // when finished with building the Ifems,
66 : // it remains to prepare the mesh for use:
67 : // find neighbors (again), partition (if needed)...
68 76 : this->_mesh.prepare_for_use ();
69 :
70 106 : return origin;
71 : }
72 :
73 :
74 :
75 :
76 :
77 :
78 :
79 :
80 :
81 :
82 :
83 :
84 10 : const Point InfElemBuilder::build_inf_elem (const InfElemOriginValue & origin_x,
85 : const InfElemOriginValue & origin_y,
86 : const InfElemOriginValue & origin_z,
87 : const bool x_sym,
88 : const bool y_sym,
89 : const bool z_sym,
90 : const bool be_verbose,
91 : std::vector<const Node *> * inner_boundary_nodes)
92 : {
93 8 : LOG_SCOPE("build_inf_elem()", "InfElemBuilder");
94 :
95 : // first determine the origin of the
96 : // infinite elements. For this, the
97 : // origin defaults to the given values,
98 : // and may be overridden when the user
99 : // provided values
100 10 : Point origin(origin_x.second, origin_y.second, origin_z.second);
101 :
102 : // when only _one_ of the origin coordinates is _not_
103 : // given, we have to determine it on our own
104 10 : if ( !origin_x.first || !origin_y.first || !origin_z.first)
105 : {
106 : // determine origin
107 0 : const BoundingBox b_box = MeshTools::create_bounding_box(_mesh);
108 0 : const Point auto_origin = (b_box.first+b_box.second)/2;
109 :
110 : // override default values, if necessary
111 0 : if (!origin_x.first)
112 0 : origin(0) = auto_origin(0);
113 : #if LIBMESH_DIM > 1
114 0 : if (!origin_y.first)
115 0 : origin(1) = auto_origin(1);
116 : #endif
117 : #if LIBMESH_DIM > 2
118 0 : if (!origin_z.first)
119 0 : origin(2) = auto_origin(2);
120 : #endif
121 :
122 0 : if (be_verbose)
123 : {
124 0 : libMesh::out << " Origin for Infinite Elements:" << std::endl;
125 :
126 0 : if (!origin_x.first)
127 0 : libMesh::out << " determined x-coordinate" << std::endl;
128 0 : if (!origin_y.first)
129 0 : libMesh::out << " determined y-coordinate" << std::endl;
130 0 : if (!origin_z.first)
131 0 : libMesh::out << " determined z-coordinate" << std::endl;
132 :
133 0 : libMesh::out << " coordinates: ";
134 0 : origin.write_unformatted(libMesh::out);
135 0 : libMesh::out << std::endl;
136 0 : }
137 0 : }
138 :
139 10 : else if (be_verbose)
140 :
141 : {
142 4 : libMesh::out << " Origin for Infinite Elements:" << std::endl;
143 4 : libMesh::out << " coordinates: ";
144 10 : origin.write_unformatted(libMesh::out);
145 4 : libMesh::out << std::endl;
146 : }
147 :
148 :
149 :
150 : // Now that we have the origin, check if the user provided an \p
151 : // inner_boundary_nodes. If so, we pass a std::set to the actual
152 : // implementation of the build_inf_elem(), so that we can convert
153 : // this to the Node * vector
154 10 : if (inner_boundary_nodes != nullptr)
155 : {
156 : // note that the std::set that we will get
157 : // from build_inf_elem() uses the index of
158 : // the element in this->_elements vector,
159 : // and the second entry is the side index
160 : // for this element. Therefore, we do _not_
161 : // need to renumber nodes and elements
162 : // prior to building the infinite elements.
163 : //
164 : // However, note that this method here uses
165 : // node id's... Do we need to renumber?
166 :
167 :
168 : // Form the list of faces of elements which finally
169 : // will tell us which nodes should receive boundary
170 : // conditions (to form the std::vector<const Node *>)
171 : std::set<std::pair<dof_id_type,
172 0 : unsigned int>> inner_faces;
173 :
174 :
175 : // build infinite elements
176 0 : this->build_inf_elem(origin,
177 : x_sym, y_sym, z_sym,
178 : be_verbose,
179 : &inner_faces);
180 :
181 0 : if (be_verbose)
182 : {
183 0 : this->_mesh.print_info();
184 0 : libMesh::out << "Data pre-processing:" << std::endl
185 0 : << " convert the <int,int> list to a Node * list..."
186 0 : << std::endl;
187 : }
188 :
189 : // First use a std::vector<dof_id_type> that holds
190 : // the global node numbers. Then sort this vector,
191 : // so that it can be made unique (no multiple occurrence
192 : // of a node), and then finally insert the Node * in
193 : // the vector inner_boundary_nodes.
194 : //
195 : // Reserve memory for the vector<> with
196 : // 4 times the size of the number of elements in the
197 : // std::set. This is a good bet for Quad4 face elements.
198 : // For higher-order elements, this probably _has_ to lead
199 : // to additional allocations...
200 : // Practice has to show how this affects performance.
201 0 : std::vector<dof_id_type> inner_boundary_node_numbers;
202 0 : inner_boundary_node_numbers.reserve(4*inner_faces.size());
203 :
204 : // Now transform the set of pairs to a list of (possibly
205 : // duplicate) global node numbers.
206 0 : for (const auto & p : inner_faces)
207 : {
208 : // build a full-ordered side element to get _all_ the base nodes
209 0 : std::unique_ptr<Elem> side(this->_mesh.elem_ref(p.first).build_side_ptr(p.second));
210 :
211 : // insert all the node numbers in inner_boundary_node_numbers
212 0 : for (const Node & node : side->node_ref_range())
213 0 : inner_boundary_node_numbers.push_back(node.id());
214 0 : }
215 :
216 :
217 : // inner_boundary_node_numbers now still holds multiple entries of
218 : // node numbers. So first sort, then unique the vector.
219 : // Note that \p std::unique only puts the new ones in
220 : // front, while to leftovers are not deleted. Instead,
221 : // it returns a pointer to the end of the unique range.
222 : //TODO:[BSK] int_ibn_size_before is not the same type as unique_size!
223 : #ifndef NDEBUG
224 0 : const std::size_t ibn_size_before = inner_boundary_node_numbers.size();
225 : #endif
226 0 : std::sort (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());
227 : auto unique_end =
228 0 : std::unique (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());
229 :
230 0 : std::size_t unique_size = std::distance(inner_boundary_node_numbers.begin(), unique_end);
231 0 : libmesh_assert_less_equal (unique_size, ibn_size_before);
232 :
233 : // Finally, create const Node * in the inner_boundary_nodes
234 : // vector. Reserve, not resize (otherwise, the push_back
235 : // would append the interesting nodes, while nullptr-nodes
236 : // live in the resize'd area...
237 0 : inner_boundary_nodes->reserve (unique_size);
238 0 : inner_boundary_nodes->clear();
239 :
240 0 : for (const auto & dof : as_range(inner_boundary_node_numbers.begin(), unique_end))
241 : {
242 0 : const Node * node = this->_mesh.node_ptr(dof);
243 0 : inner_boundary_nodes->push_back(node);
244 : }
245 :
246 0 : if (be_verbose)
247 0 : libMesh::out << " finished identifying " << unique_size
248 0 : << " target nodes." << std::endl;
249 : }
250 :
251 : else
252 :
253 : {
254 : // There are no inner boundary nodes, so simply build the infinite elements
255 10 : this->build_inf_elem(origin, x_sym, y_sym, z_sym, be_verbose);
256 : }
257 :
258 : // when finished with building the Ifems,
259 : // it remains to prepare the mesh for use:
260 : // find neighbors again, partition (if needed)...
261 10 : this->_mesh.prepare_for_use ();
262 :
263 14 : return origin;
264 : }
265 :
266 :
267 :
268 :
269 :
270 :
271 :
272 :
273 :
274 : // The actual implementation of building elements.
275 86 : void InfElemBuilder::build_inf_elem(const Point & origin,
276 : const bool x_sym,
277 : const bool y_sym,
278 : const bool z_sym,
279 : const bool be_verbose,
280 : std::set<std::pair<dof_id_type,
281 : unsigned int>> * inner_faces)
282 : {
283 86 : if (be_verbose)
284 : {
285 : #ifdef DEBUG
286 8 : libMesh::out << " Building Infinite Elements:" << std::endl;
287 8 : libMesh::out << " updating element neighbor tables..." << std::endl;
288 : #else
289 : libMesh::out << " Verbose mode disabled in non-debug mode." << std::endl;
290 : #endif
291 : }
292 :
293 :
294 : // update element neighbors
295 86 : this->_mesh.find_neighbors();
296 :
297 68 : LOG_SCOPE("build_inf_elem()", "InfElemBuilder");
298 :
299 : // A set for storing element number, side number pairs.
300 : // pair.first == element number, pair.second == side number
301 68 : std::set<std::pair<dof_id_type,unsigned int>> faces;
302 68 : std::set<std::pair<dof_id_type,unsigned int>> ofaces;
303 :
304 : // A set for storing node numbers on the outer faces.
305 68 : std::set<dof_id_type> onodes;
306 :
307 : // The distance to the farthest point in the mesh from the origin
308 86 : Real max_r=0.;
309 :
310 : // The index of the farthest point in the mesh from the origin
311 86 : int max_r_node = -1;
312 :
313 : #ifdef DEBUG
314 34 : if (be_verbose)
315 : {
316 8 : libMesh::out << " collecting boundary sides";
317 8 : if (x_sym || y_sym || z_sym)
318 0 : libMesh::out << ", skipping sides in symmetry planes..." << std::endl;
319 : else
320 8 : libMesh::out << "..." << std::endl;
321 : }
322 : #endif
323 :
324 : // Iterate through all elements and sides, collect indices of all active
325 : // boundary sides in the faces set. Skip sides which lie in symmetry planes.
326 : // Later, sides of the inner boundary will be sorted out.
327 8600 : for (const auto & elem : _mesh.active_element_ptr_range())
328 44147 : for (auto s : elem->side_index_range())
329 47984 : if (elem->neighbor_ptr(s) == nullptr)
330 : {
331 : // note that it is safe to use the Elem::side() method,
332 : // which gives a non-full-ordered element
333 4934 : std::unique_ptr<Elem> side(elem->build_side_ptr(s));
334 :
335 : // bool flags for symmetry detection
336 1408 : bool sym_side=false;
337 1408 : bool on_x_sym=true;
338 1408 : bool on_y_sym=true;
339 1408 : bool on_z_sym=true;
340 :
341 :
342 : // Loop over the nodes to check whether they are on the symmetry planes,
343 : // and therefore sufficient to use a non-full-ordered side element
344 25570 : for (const Node & node : side->node_ref_range())
345 : {
346 17616 : const dof_id_type node_id = node.id();
347 : const Point dist_from_origin =
348 22044 : this->_mesh.point(node_id) - origin;
349 :
350 22044 : if (x_sym)
351 0 : if (std::abs(dist_from_origin(0)) > 1.e-3)
352 0 : on_x_sym=false;
353 :
354 22044 : if (y_sym)
355 0 : if (std::abs(dist_from_origin(1)) > 1.e-3)
356 0 : on_y_sym=false;
357 :
358 22044 : if (z_sym)
359 0 : if (std::abs(dist_from_origin(2)) > 1.e-3)
360 0 : on_z_sym=false;
361 :
362 : // if (x_sym)
363 : // if (std::abs(dist_from_origin(0)) > 1.e-6)
364 : // on_x_sym=false;
365 :
366 : // if (y_sym)
367 : // if (std::abs(dist_from_origin(1)) > 1.e-6)
368 : // on_y_sym=false;
369 :
370 : // if (z_sym)
371 : // if (std::abs(dist_from_origin(2)) > 1.e-6)
372 : // on_z_sym=false;
373 :
374 : //find the node most distant from origin
375 :
376 13236 : Real r = dist_from_origin.norm();
377 22044 : if (r > max_r)
378 : {
379 157 : max_r = r;
380 157 : max_r_node=node_id;
381 : }
382 :
383 : }
384 :
385 3526 : sym_side = (x_sym && on_x_sym) || (y_sym && on_y_sym) || (z_sym && on_z_sym);
386 :
387 1408 : if (!sym_side)
388 3526 : faces.emplace(elem->id(), s);
389 :
390 728 : } // neighbor(s) == nullptr
391 :
392 : // On a distributed mesh it might be some other processor who sees
393 : // the farthest node.
394 86 : if (!this->_mesh.is_serial())
395 : {
396 : unsigned int rank;
397 0 : this->_mesh.comm().maxloc(max_r, rank);
398 0 : this->_mesh.comm().broadcast(max_r_node, rank);
399 : }
400 :
401 : // If a boundary side has one node on the outer boundary,
402 : // all points of this side are on the outer boundary.
403 : // Start with the node most distant from origin, which has
404 : // to be on the outer boundary, then recursively find all
405 : // sides and nodes connected to it. Found sides are moved
406 : // from faces to ofaces, nodes are collected in onodes.
407 : // Here, the search is done iteratively, because, depending on
408 : // the mesh, a very high level of recursion might be necessary.
409 86 : if (max_r_node >= 0)
410 : // include the possibility of the 1st element being most far away.
411 : // Only the case of no outer boundary is to be excluded.
412 86 : onodes.insert(max_r_node);
413 :
414 : // If we're not on a serial mesh, we'll need to synchronize that
415 : // onodes list too.
416 : bool did_parallel_update;
417 :
418 86 : do
419 : {
420 86 : did_parallel_update = false;
421 :
422 34 : auto face_it = faces.begin();
423 34 : auto face_end = faces.end();
424 34 : unsigned int facesfound=0;
425 4712 : while (face_it != face_end) {
426 4626 : std::pair<dof_id_type, unsigned int> p = *face_it;
427 :
428 : // This has to be a full-ordered side element,
429 : // since we need the correct n_nodes,
430 6474 : std::unique_ptr<Elem> side(this->_mesh.elem_ref(p.first).build_side_ptr(p.second));
431 :
432 1848 : bool found=false;
433 13127 : for (const Node & node : side->node_ref_range())
434 12027 : if (onodes.count(node.id()))
435 : {
436 1408 : found=true;
437 1408 : break;
438 : }
439 :
440 : // If a new oface is found, include its nodes in onodes
441 4626 : if (found)
442 : {
443 25570 : for (const Node & node : side->node_ref_range())
444 22044 : onodes.insert(node.id());
445 :
446 2118 : ofaces.insert(p);
447 2118 : face_it = faces.erase(face_it); // increment is done here
448 :
449 3526 : facesfound++;
450 : }
451 :
452 : else
453 440 : ++face_it; // increment is done here
454 :
455 : // If at least one new oface was found in this cycle,
456 : // do another search cycle.
457 4626 : if (facesfound>0 && face_it == faces.end())
458 : {
459 44 : facesfound = 0;
460 44 : face_it = faces.begin();
461 : }
462 930 : }
463 :
464 86 : if (!this->_mesh.is_serial())
465 : {
466 0 : auto my_onodes_size = onodes.size();
467 0 : this->_mesh.comm().set_union(onodes);
468 0 : did_parallel_update = (onodes.size() > my_onodes_size);
469 0 : this->_mesh.comm().max(did_parallel_update);
470 : }
471 : }
472 : while (did_parallel_update);
473 :
474 :
475 : #ifdef DEBUG
476 34 : if (be_verbose)
477 8 : libMesh::out << " found "
478 8 : << faces.size()
479 8 : << " inner and "
480 16 : << ofaces.size()
481 8 : << " outer boundary faces"
482 8 : << std::endl;
483 : #endif
484 :
485 : // When the user provided a non-null pointer to
486 : // inner_faces, that implies he wants to have
487 : // this std::set. For now, simply copy the data.
488 86 : if (inner_faces != nullptr)
489 0 : *inner_faces = faces;
490 :
491 : // free memory, clear our local variable, no need
492 : // for it any more.
493 34 : faces.clear();
494 :
495 :
496 : // outer_nodes maps onodes to their duplicates
497 68 : std::map<dof_id_type, Node *> outer_nodes;
498 :
499 : // We may need to pick our own object ids in parallel
500 86 : dof_id_type old_max_node_id = _mesh.max_node_id();
501 86 : dof_id_type old_max_elem_id = _mesh.max_elem_id();
502 :
503 : // Likewise with our unique_ids
504 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
505 86 : unique_id_type old_max_unique_id = _mesh.parallel_max_unique_id();
506 : #endif
507 :
508 : // for each boundary node, add an outer_node with
509 : // double distance from origin.
510 8184 : for (const auto & n : onodes)
511 : {
512 8098 : if (!this->_mesh.query_node_ptr(n))
513 : {
514 0 : libmesh_assert(!_mesh.is_serial());
515 0 : continue;
516 : }
517 :
518 8098 : Point p = (Point(this->_mesh.point(n)) * 2) - origin;
519 8098 : if (_mesh.is_serial())
520 : {
521 : // Add with a default id in serial
522 8098 : outer_nodes[n]=this->_mesh.add_point(p);
523 : }
524 : else
525 : {
526 : // Pick a unique id in parallel
527 0 : Node & bnode = _mesh.node_ref(n);
528 0 : dof_id_type new_id = bnode.id() + old_max_node_id;
529 0 : std::unique_ptr<Node> new_node = Node::build(p, new_id);
530 0 : new_node->processor_id() = bnode.processor_id();
531 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
532 0 : new_node->set_unique_id(old_max_unique_id + bnode.id());
533 : #endif
534 :
535 0 : outer_nodes[n] =
536 0 : this->_mesh.add_node(std::move(new_node));
537 0 : }
538 : }
539 :
540 :
541 : #ifdef DEBUG
542 : // for verbose, remember n_elem
543 34 : dof_id_type n_conventional_elem = this->_mesh.n_elem();
544 : #endif
545 :
546 :
547 : // build Elems based on boundary side type
548 3612 : for (auto & p : ofaces)
549 : {
550 3526 : Elem & belem = this->_mesh.elem_ref(p.first);
551 :
552 : // build a full-ordered side element to get the base nodes
553 3526 : std::unique_ptr<Elem> side(belem.build_side_ptr(p.second));
554 :
555 : // create cell depending on side type, assign nodes,
556 : // use braces to force scope.
557 1408 : bool is_higher_order_elem = false;
558 :
559 710 : std::unique_ptr<Elem> el;
560 3526 : switch(side->type())
561 : {
562 : // 3D infinite elements
563 : // TRIs
564 0 : case TRI3:
565 0 : el = Elem::build(INFPRISM6);
566 0 : break;
567 :
568 1820 : case TRI6:
569 2184 : el = Elem::build(INFPRISM12);
570 728 : is_higher_order_elem = true;
571 1820 : break;
572 :
573 : // QUADs
574 846 : case QUAD4:
575 1020 : el = Elem::build(INFHEX8);
576 846 : break;
577 :
578 0 : case QUAD8:
579 0 : el = Elem::build(INFHEX16);
580 0 : is_higher_order_elem = true;
581 0 : break;
582 :
583 860 : case QUAD9:
584 1032 : el = Elem::build(INFHEX18);
585 :
586 : // the method of assigning nodes (which follows below)
587 : // omits in the case of QUAD9 the bubble node; therefore
588 : // we assign these first by hand here.
589 1548 : el->set_node(16, side->node_ptr(8));
590 1548 : el->set_node(17, outer_nodes[side->node_id(8)]);
591 344 : is_higher_order_elem=true;
592 860 : break;
593 :
594 : // 2D infinite elements
595 0 : case EDGE2:
596 0 : el = Elem::build(INFQUAD4);
597 0 : break;
598 :
599 0 : case EDGE3:
600 0 : el = Elem::build(INFQUAD6);
601 0 : el->set_node(4, side->node_ptr(2));
602 0 : break;
603 :
604 : // 1D infinite elements not supported
605 0 : default:
606 0 : libMesh::out << "InfElemBuilder::build_inf_elem(Point, bool, bool, bool, bool): "
607 0 : << "invalid face element "
608 0 : << std::endl;
609 0 : continue;
610 : }
611 :
612 3526 : const unsigned int n_base_vertices = side->n_vertices();
613 :
614 : // On a distributed mesh, manually assign unique ids to the new
615 : // element, and make sure any RemoteElem neighbor links are set.
616 3526 : if (!_mesh.is_serial())
617 : {
618 0 : el->processor_id() = belem.processor_id();
619 :
620 : // We'd better not have elements with more than 6 sides
621 0 : const unsigned int max_sides = 6;
622 0 : libmesh_assert_less_equal(el->n_sides(), max_sides);
623 0 : el->set_id (belem.id() * max_sides + p.second + old_max_elem_id);
624 :
625 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
626 0 : el->set_unique_id(old_max_unique_id + old_max_node_id +
627 0 : belem.id() * max_sides + p.second);
628 : #endif
629 :
630 : // If we have a remote neighbor on a boundary element side
631 0 : if (belem.dim() > 1)
632 0 : for (auto s : belem.side_index_range())
633 0 : if (belem.neighbor_ptr(s) == remote_elem)
634 : {
635 : // Find any corresponding infinite element side
636 0 : std::unique_ptr<const Elem> remote_side(belem.build_side_ptr(s));
637 :
638 0 : for (auto inf_s : el->side_index_range())
639 : {
640 : // The base side 0 shares all vertices but isn't
641 : // remote
642 0 : if (!inf_s)
643 0 : continue;
644 :
645 : // But another side, one which shares enough
646 : // vertices to show it's the same side, is.
647 0 : unsigned int n_shared_vertices = 0;
648 0 : for (unsigned int i = 0; i != n_base_vertices; ++i)
649 0 : for (auto & node : remote_side->node_ref_range())
650 0 : if (side->node_ptr(i) == &node &&
651 0 : el->is_node_on_side(i,inf_s))
652 0 : ++n_shared_vertices;
653 :
654 0 : if (n_shared_vertices + 1 >= belem.dim())
655 : {
656 : el->set_neighbor
657 0 : (inf_s, const_cast<RemoteElem *>(remote_elem));
658 0 : break;
659 : }
660 : }
661 0 : }
662 : }
663 :
664 : // assign vertices to the new infinite element
665 15810 : for (unsigned int i=0; i<n_base_vertices; i++)
666 : {
667 22092 : el->set_node(i , side->node_ptr(i));
668 22092 : el->set_node(i+n_base_vertices, outer_nodes[side->node_id(i)]);
669 : }
670 :
671 :
672 : // when this is a higher order element,
673 : // assign also the nodes in between
674 3526 : if (is_higher_order_elem)
675 : {
676 : // n_safe_base_nodes is the number of nodes in \p side
677 : // that may be safely assigned using below for loop.
678 : // Actually, n_safe_base_nodes is _identical_ with el->n_vertices(),
679 : // since for QUAD9, the 9th node was already assigned above
680 2680 : const unsigned int n_safe_base_nodes = el->n_vertices();
681 :
682 11580 : for (unsigned int i=n_base_vertices; i<n_safe_base_nodes; i++)
683 : {
684 16020 : el->set_node(i+n_base_vertices, side->node_ptr(i));
685 8900 : el->set_node(i+n_safe_base_nodes,
686 21360 : outer_nodes[side->node_id(i)]);
687 : }
688 : }
689 :
690 :
691 : // add infinite element to mesh
692 6342 : this->_mesh.add_elem(std::move(el));
693 710 : } // for
694 :
695 : // If we're not in serial, we might not have all our remote_elem
696 : // neighbors set up properly on the new elements: two new infinite
697 : // elements which should be neighbors may not be rooted in finite
698 : // elmeents which are neighbors, and we were only looking for
699 : // remote_elem links on each root finite element. We'll fix the
700 : // problem by first finding all the neighbors we can and then
701 : // recognizing that any missing neighbors on infinite elements must
702 : // be remote.
703 86 : if (!_mesh.is_serial())
704 : {
705 0 : _mesh.find_neighbors(/*reset_remote_elements=*/ false,
706 : /*reset_current_list=*/ true,
707 0 : /*assert_valid=*/ false);
708 :
709 0 : for (auto & p : ofaces)
710 : {
711 0 : Elem & belem = this->_mesh.elem_ref(p.first);
712 0 : Elem * inf_elem = belem.neighbor_ptr(p.second);
713 0 : libmesh_assert(inf_elem);
714 :
715 0 : for (auto s : make_range(inf_elem->n_sides()))
716 0 : if (!inf_elem->neighbor_ptr(s))
717 : inf_elem->set_neighbor
718 0 : (s, const_cast<RemoteElem *>(remote_elem));
719 : }
720 : }
721 :
722 : #ifdef DEBUG
723 34 : _mesh.libmesh_assert_valid_parallel_ids();
724 :
725 34 : if (be_verbose)
726 8 : libMesh::out << " added "
727 8 : << this->_mesh.n_elem() - n_conventional_elem
728 8 : << " infinite elements and "
729 16 : << onodes.size()
730 8 : << " nodes to the mesh"
731 8 : << std::endl
732 8 : << std::endl;
733 : #endif
734 86 : }
735 :
736 : } // namespace libMesh
737 :
738 :
739 :
740 :
741 :
742 : #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
|
