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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo