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