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