serial_mesh.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 
19 
20 // Local includes
21 #include "libmesh/boundary_info.h"
22 #include "libmesh/elem.h"
25 #include "libmesh/serial_mesh.h"
26 #include "libmesh/utility.h"
27 
28 #include LIBMESH_INCLUDE_UNORDERED_SET
29 LIBMESH_DEFINE_HASH_POINTERS
30 
31 
32 namespace
33 {
34  // A custom comparison function, based on Point::operator<,
35  // that tries to ignore floating point differences in components
36  // of the point
38  {
39  private:
41 
42  public:
43  // Constructor takes the tolerance to be used in fuzzy comparisons
44  FuzzyPointCompare(Real tol) : _tol(tol) {}
45 
46  // This is inspired directly by Point::operator<
47  bool operator()(const Point& lhs, const Point& rhs)
48  {
49  for (unsigned i=0; i<LIBMESH_DIM; ++i)
50  {
51  // If the current components are within some tolerance
52  // of one another, then don't attempt the less-than comparison.
53  // Note that this may cause something strange to happen, as Roy
54  // believes he can prove it is not a total ordering...
55  Real rel_size = std::max(std::abs(lhs(i)), std::abs(rhs(i)));
56 
57  // Don't use relative tolerance if both numbers are already small.
58  // How small? Some possible options are:
59  // * std::numeric_limits<Real>::epsilon()
60  // * TOLERANCE
61  // * 1.0
62  // If we use std::numeric_limits<Real>::epsilon(), we'll
63  // do more relative comparisons for small numbers, but
64  // increase the chance for false positives? If we pick 1.0,
65  // we'll never "increase" the difference between small numbers
66  // in the test below.
67  if (rel_size < 1.)
68  rel_size = 1.;
69 
70  // Don't attempt the comparison if lhs(i) and rhs(i) are too close
71  // together.
72  if ( std::abs(lhs(i) - rhs(i)) / rel_size < _tol)
73  continue;
74 
75  if (lhs(i) < rhs(i))
76  return true;
77  if (lhs(i) > rhs(i))
78  return false;
79  }
80 
81  // We compared all the components without returning yet, so
82  // each component was neither greater than nor less than they other.
83  // They might be equal, so return false.
84  return false;
85  }
86 
87  // Needed by std::sort on vector< pair<Point,id> >
88  bool operator()(const std::pair<Point, dof_id_type>& lhs,
89  const std::pair<Point, dof_id_type>& rhs)
90  {
91  return (*this)(lhs.first, rhs.first);
92  }
93 
94  // Comparsion function where lhs is a Point and rhs is a pair<Point,dof_id_type>.
95  // This is used in routines like lower_bound, where a specific value is being
96  // searched for.
97  bool operator()(const Point& lhs, std::pair<Point, dof_id_type>& rhs)
98  {
99  return (*this)(lhs, rhs.first);
100  }
101 
102  // And the other way around...
103  bool operator()(std::pair<Point, dof_id_type>& lhs, const Point& rhs)
104  {
105  return (*this)(lhs.first, rhs);
106  }
107  };
108 }
109 
110 
111 
112 namespace libMesh
113 {
114 
115 // ------------------------------------------------------------
116 // SerialMesh class member functions
118  unsigned int d) :
119  UnstructuredMesh (comm,d)
120 {
121 #ifdef LIBMESH_ENABLE_UNIQUE_ID
122  // In serial we just need to reset the next unique id to zero
123  // here in the constructor.
124  _next_unique_id = 0;
125 #endif
127 }
128 
129 
130 
131 #ifndef LIBMESH_DISABLE_COMMWORLD
132 SerialMesh::SerialMesh (unsigned int d) :
133  UnstructuredMesh (d)
134 {
135 #ifdef LIBMESH_ENABLE_UNIQUE_ID
136  // In serial we just need to reset the next unique id to zero
137  // here in the constructor.
138  _next_unique_id = 0;
139 #endif
141 }
142 #endif
143 
144 
146 {
147  this->clear(); // Free nodes and elements
148 }
149 
150 
151 // This might be specialized later, but right now it's just here to
152 // make sure the compiler doesn't give us a default (non-deep) copy
153 // constructor instead.
154 SerialMesh::SerialMesh (const SerialMesh &other_mesh) :
155  UnstructuredMesh (other_mesh)
156 {
157  this->copy_nodes_and_elements(other_mesh);
158  *this->boundary_info = *other_mesh.boundary_info;
159 }
160 
161 
163  UnstructuredMesh (other_mesh)
164 {
165  this->copy_nodes_and_elements(other_mesh);
166  *this->boundary_info = *other_mesh.boundary_info;
167 }
168 
169 
170 const Point& SerialMesh::point (const dof_id_type i) const
171 {
172  libmesh_assert_less (i, this->n_nodes());
174  libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
175 
176  return (*_nodes[i]);
177 }
178 
179 
180 
181 
182 
183 const Node& SerialMesh::node (const dof_id_type i) const
184 {
185  libmesh_assert_less (i, this->n_nodes());
187  libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
188 
189  return (*_nodes[i]);
190 }
191 
192 
193 
194 
195 
197 {
198  if (i >= this->n_nodes())
199  {
200  libMesh::out << " i=" << i
201  << ", n_nodes()=" << this->n_nodes()
202  << std::endl;
203  libmesh_error();
204  }
205 
206  libmesh_assert_less (i, this->n_nodes());
208  libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
209 
210  return (*_nodes[i]);
211 }
212 
213 
214 
215 const Node* SerialMesh::node_ptr (const dof_id_type i) const
216 {
217  libmesh_assert_less (i, this->n_nodes());
219  libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
220 
221  return _nodes[i];
222 }
223 
224 
225 
226 
228 {
229  libmesh_assert_less (i, this->n_nodes());
231  libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
232 
233  return _nodes[i];
234 }
235 
236 
237 
238 
240 {
241  if (i >= this->n_nodes())
242  return NULL;
243  libmesh_assert (_nodes[i] == NULL ||
244  _nodes[i]->id() == i); // This will change soon
245 
246  return _nodes[i];
247 }
248 
249 
250 
251 
253 {
254  if (i >= this->n_nodes())
255  return NULL;
256  libmesh_assert (_nodes[i] == NULL ||
257  _nodes[i]->id() == i); // This will change soon
258 
259  return _nodes[i];
260 }
261 
262 
263 
264 
265 const Elem* SerialMesh::elem (const dof_id_type i) const
266 {
267  libmesh_assert_less (i, this->n_elem());
269  libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
270 
271  return _elements[i];
272 }
273 
274 
275 
276 
278 {
279  libmesh_assert_less (i, this->n_elem());
281  libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
282 
283  return _elements[i];
284 }
285 
286 
287 
288 
290 {
291  if (i >= this->n_elem())
292  return NULL;
293  libmesh_assert (_elements[i] == NULL ||
294  _elements[i]->id() == i); // This will change soon
295 
296  return _elements[i];
297 }
298 
299 
300 
301 
303 {
304  if (i >= this->n_elem())
305  return NULL;
306  libmesh_assert (_elements[i] == NULL ||
307  _elements[i]->id() == i); // This will change soon
308 
309  return _elements[i];
310 }
311 
312 
313 
314 
316 {
317  libmesh_assert(e);
318 
319  // We no longer merely append elements with SerialMesh
320 
321  // If the user requests a valid id that doesn't correspond to an
322  // existing element, let's give them that id, resizing the elements
323  // container if necessary.
324  if (!e->valid_id())
325  e->set_id (_elements.size());
326 
327 #ifdef LIBMESH_ENABLE_UNIQUE_ID
328  if (!e->valid_unique_id())
330 #endif
331 
332  const dof_id_type id = e->id();
333 
334  if (id < _elements.size())
335  {
336  // Overwriting existing elements is still probably a mistake.
338  }
339  else
340  {
341  _elements.resize(id+1, NULL);
342  }
343 
344  _elements[id] = e;
345 
346  return e;
347 }
348 
349 
350 
352 {
353  dof_id_type eid = e->id();
354  libmesh_assert_less (eid, _elements.size());
355  Elem *oldelem = _elements[eid];
356 
357  if (oldelem)
358  {
359  libmesh_assert_equal_to (oldelem->id(), eid);
360  this->delete_elem(oldelem);
361  }
362 
363  _elements[e->id()] = e;
364 
365  return e;
366 }
367 
368 
369 
371 {
372  libmesh_assert(e);
373 
374  // Initialize an iterator to eventually point to the element we want to delete
375  std::vector<Elem*>::iterator pos = _elements.end();
376 
377  // In many cases, e->id() gives us a clue as to where e
378  // is located in the _elements vector. Try that first
379  // before trying the O(n_elem) search.
380  libmesh_assert_less (e->id(), _elements.size());
381 
382  if (_elements[e->id()] == e)
383  {
384  // We found it!
385  pos = _elements.begin();
386  std::advance(pos, e->id());
387  }
388 
389  else
390  {
391  // This search is O(n_elem)
392  pos = std::find (_elements.begin(),
393  _elements.end(),
394  e);
395  }
396 
397  // Huh? Element not in the vector?
398  libmesh_assert (pos != _elements.end());
399 
400  // Remove the element from the BoundaryInfo object
401  this->boundary_info->remove(e);
402 
403  // delete the element
404  delete e;
405 
406  // explicitly NULL the pointer
407  *pos = NULL;
408 }
409 
410 
411 
413  const dof_id_type new_id)
414 {
415  // This doesn't get used in serial yet
416  Elem *el = _elements[old_id];
417  libmesh_assert (el);
418 
419  el->set_id(new_id);
420  libmesh_assert (!_elements[new_id]);
421  _elements[new_id] = el;
422  _elements[old_id] = NULL;
423 }
424 
425 
426 
428  const dof_id_type id,
429  const processor_id_type proc_id)
430 {
431 // // We only append points with SerialMesh
432 // libmesh_assert(id == DofObject::invalid_id || id == _nodes.size());
433 // Node *n = Node::build(p, _nodes.size()).release();
434 // n->processor_id() = proc_id;
435 // _nodes.push_back (n);
436 
437  Node *n = NULL;
438 
439  // If the user requests a valid id, either
440  // provide the existing node or resize the container
441  // to fit the new node.
442  if (id != DofObject::invalid_id)
443  if (id < _nodes.size())
444  n = _nodes[id];
445  else
446  _nodes.resize(id+1);
447  else
448  _nodes.push_back (static_cast<Node*>(NULL));
449 
450  // if the node already exists, then assign new (x,y,z) values
451  if (n)
452  *n = p;
453  // otherwise build a new node, put it in the right spot, and return
454  // a valid pointer.
455  else
456  {
457  n = Node::build(p, (id == DofObject::invalid_id) ? _nodes.size()-1 : id).release();
458  n->processor_id() = proc_id;
459 
460  if (id == DofObject::invalid_id)
461  _nodes.back() = n;
462  else
463  _nodes[id] = n;
464  }
465 
466  // better not pass back a NULL pointer.
467  libmesh_assert (n);
468 
469  return n;
470 }
471 
472 
473 
475 {
476  libmesh_assert(n);
477  // We only append points with SerialMesh
478  libmesh_assert(!n->valid_id() || n->id() == _nodes.size());
479 
480  n->set_id (_nodes.size());
481 
482 #ifdef LIBMESH_ENABLE_UNIQUE_ID
483  if (!n->valid_unique_id())
485 #endif
486 
487  _nodes.push_back(n);
488 
489  return n;
490 }
491 
492 
493 
495 {
496  if (!n)
497  {
498  libMesh::err << "Error, attempting to insert NULL node." << std::endl;
499  libmesh_error();
500  }
501 
502  if (n->id() == DofObject::invalid_id)
503  {
504  libMesh::err << "Error, cannot insert node with invalid id." << std::endl;
505  libmesh_error();
506  }
507 
508  if (n->id() < _nodes.size())
509  {
510  // Don't allow inserting on top of an existing Node.
511  if (_nodes[ n->id() ] != NULL)
512  {
513  libMesh::err << "Error, cannot insert node on top of existing node." << std::endl;
514  libmesh_error();
515  }
516  }
517  else
518  {
519  // Allocate just enough space to store the new node. This will
520  // cause highly non-ideal memory allocation behavior if called
521  // repeatedly...
522  _nodes.resize(n->id() + 1);
523  }
524 
525 
526  // We have enough space and this spot isn't already occupied by
527  // another node, so go ahead and add it.
528  _nodes[ n->id() ] = n;
529 
530  // If we made it this far, we just inserted the node the user handed
531  // us, so we can give it right back.
532  return n;
533 }
534 
535 
536 
538 {
539  libmesh_assert(n);
540  libmesh_assert_less (n->id(), _nodes.size());
541 
542  // Initialize an iterator to eventually point to the element we want
543  // to delete
544  std::vector<Node*>::iterator pos;
545 
546  // In many cases, e->id() gives us a clue as to where e
547  // is located in the _elements vector. Try that first
548  // before trying the O(n_elem) search.
549  if (_nodes[n->id()] == n)
550  {
551  pos = _nodes.begin();
552  std::advance(pos, n->id());
553  }
554  else
555  {
556  pos = std::find (_nodes.begin(),
557  _nodes.end(),
558  n);
559  }
560 
561  // Huh? Node not in the vector?
562  libmesh_assert (pos != _nodes.end());
563 
564  // Delete the node from the BoundaryInfo object
565  this->boundary_info->remove(n);
566 
567  // delete the node
568  delete n;
569 
570  // explicitly NULL the pointer
571  *pos = NULL;
572 }
573 
574 
575 
577  const dof_id_type new_id)
578 {
579  // This doesn't get used in serial yet
580  Node *nd = _nodes[old_id];
581  libmesh_assert (nd);
582 
583  nd->set_id(new_id);
584  libmesh_assert (!_nodes[new_id]);
585  _nodes[new_id] = nd;
586  _nodes[old_id] = NULL;
587 }
588 
589 
590 
592 {
593  // Call parent clear function
594  MeshBase::clear();
595 
596 
597  // Clear our elements and nodes
598  {
599  std::vector<Elem*>::iterator it = _elements.begin();
600  const std::vector<Elem*>::iterator end = _elements.end();
601 
602  // There is no need to remove the elements from
603  // the BoundaryInfo data structure since we
604  // already cleared it.
605  for (; it != end; ++it)
606  delete *it;
607 
608  _elements.clear();
609  }
610 
611  // clear the nodes data structure
612  {
613  std::vector<Node*>::iterator it = _nodes.begin();
614  const std::vector<Node*>::iterator end = _nodes.end();
615 
616  // There is no need to remove the nodes from
617  // the BoundaryInfo data structure since we
618  // already cleared it.
619  for (; it != end; ++it)
620  delete *it;
621 
622  _nodes.clear();
623  }
624 }
625 
626 
627 
629 {
630 
631  START_LOG("renumber_nodes_and_elem()", "Mesh");
632 
633  // node and element id counters
634  dof_id_type next_free_elem = 0;
635  dof_id_type next_free_node = 0;
636 
637  // Will hold the set of nodes that are currently connected to elements
638  LIBMESH_BEST_UNORDERED_SET<Node*> connected_nodes;
639 
640  // Loop over the elements. Note that there may
641  // be NULLs in the _elements vector from the coarsening
642  // process. Pack the elements in to a contiguous array
643  // and then trim any excess.
644  {
645  std::vector<Elem*>::iterator in = _elements.begin();
646  std::vector<Elem*>::iterator out_iter = _elements.begin();
647  const std::vector<Elem*>::iterator end = _elements.end();
648 
649  for (; in != end; ++in)
650  if (*in != NULL)
651  {
652  Elem* el = *in;
653 
654  *out_iter = *in;
655  ++out_iter;
656 
657  // Increment the element counter
658  el->set_id (next_free_elem++);
659 
661  {
662  // Add this elements nodes to the connected list
663  for (unsigned int n=0; n<el->n_nodes(); n++)
664  connected_nodes.insert(el->get_node(n));
665  }
666  else // We DO want node renumbering
667  {
668  // Loop over this element's nodes. Number them,
669  // if they have not been numbered already. Also,
670  // position them in the _nodes vector so that they
671  // are packed contiguously from the beginning.
672  for (unsigned int n=0; n<el->n_nodes(); n++)
673  if (el->node(n) == next_free_node) // don't need to process
674  next_free_node++; // [(src == dst) below]
675 
676  else if (el->node(n) > next_free_node) // need to process
677  {
678  // The source and destination indices
679  // for this node
680  const dof_id_type src_idx = el->node(n);
681  const dof_id_type dst_idx = next_free_node++;
682 
683  // ensure we want to swap a valid nodes
684  libmesh_assert(_nodes[src_idx]);
685 
686  // Swap the source and destination nodes
687  std::swap(_nodes[src_idx],
688  _nodes[dst_idx] );
689 
690  // Set proper indices where that makes sense
691  if (_nodes[src_idx] != NULL)
692  _nodes[src_idx]->set_id (src_idx);
693  _nodes[dst_idx]->set_id (dst_idx);
694  }
695  }
696  }
697 
698  // Erase any additional storage. These elements have been
699  // copied into NULL voids by the procedure above, and are
700  // thus repeated and unnecessary.
701  _elements.erase (out_iter, end);
702  }
703 
704 
706  {
707  // Loop over the nodes. Note that there may
708  // be NULLs in the _nodes vector from the coarsening
709  // process. Pack the nodes in to a contiguous array
710  // and then trim any excess.
711 
712  std::vector<Node*>::iterator in = _nodes.begin();
713  std::vector<Node*>::iterator out_iter = _nodes.begin();
714  const std::vector<Node*>::iterator end = _nodes.end();
715 
716  for (; in != end; ++in)
717  if (*in != NULL)
718  {
719  // This is a reference so that if we change the pointer it will change in the vector
720  Node* & nd = *in;
721 
722  // If this node is still connected to an elem, put it in the list
723  if(connected_nodes.find(nd) != connected_nodes.end())
724  {
725  *out_iter = nd;
726  ++out_iter;
727 
728  // Increment the node counter
729  nd->set_id (next_free_node++);
730  }
731  else // This node is orphaned, delete it!
732  {
733  this->boundary_info->remove (nd);
734 
735  // delete the node
736  delete nd;
737  nd = NULL;
738  }
739  }
740 
741  // Erase any additional storage. Whatever was
742  _nodes.erase (out_iter, end);
743  }
744  else // We really DO want node renumbering
745  {
746  // Any nodes in the vector >= _nodes[next_free_node]
747  // are not connected to any elements and may be deleted
748  // if desired.
749 
750  // (This code block will erase the unused nodes)
751  // Now, delete the unused nodes
752  {
753  std::vector<Node*>::iterator nd = _nodes.begin();
754  const std::vector<Node*>::iterator end = _nodes.end();
755 
756  std::advance (nd, next_free_node);
757 
758  for (std::vector<Node*>::iterator it=nd;
759  it != end; ++it)
760  {
761  // Mesh modification code might have already deleted some
762  // nodes
763  if (*it == NULL)
764  continue;
765 
766  // remove any boundary information associated with
767  // this node
768  this->boundary_info->remove (*it);
769 
770  // delete the node
771  delete *it;
772  *it = NULL;
773  }
774 
775  _nodes.erase (nd, end);
776  }
777  }
778 
779  libmesh_assert_equal_to (next_free_elem, _elements.size());
780  libmesh_assert_equal_to (next_free_node, _nodes.size());
781 
782  STOP_LOG("renumber_nodes_and_elem()", "Mesh");
783 }
784 
785 
786 
788 {
789  // Nodes first
790  for (dof_id_type n=0; n<this->_nodes.size(); n++)
791  if (this->_nodes[n] != NULL)
792  this->_nodes[n]->set_id() = n;
793 
794  // Elements next
795  for (dof_id_type e=0; e<this->_elements.size(); e++)
796  if (this->_elements[e] != NULL)
797  this->_elements[e]->set_id() = e;
798 }
799 
800 
802  boundary_id_type this_mesh_boundary_id,
803  boundary_id_type other_mesh_boundary_id,
804  Real tol,
805  bool clear_stitched_boundary_ids,
806  bool verbose,
807  bool use_binary_search,
808  bool enforce_all_nodes_match_on_boundaries)
809 {
810  stitching_helper(&other_mesh,
811  this_mesh_boundary_id,
812  other_mesh_boundary_id,
813  tol,
814  clear_stitched_boundary_ids,
815  verbose,
816  use_binary_search,
817  enforce_all_nodes_match_on_boundaries);
818 }
819 
821  boundary_id_type boundary_id_2,
822  Real tol,
823  bool clear_stitched_boundary_ids,
824  bool verbose,
825  bool use_binary_search,
826  bool enforce_all_nodes_match_on_boundaries)
827 {
828  stitching_helper(NULL,
829  boundary_id_1,
830  boundary_id_2,
831  tol,
832  clear_stitched_boundary_ids,
833  verbose,
834  use_binary_search,
835  enforce_all_nodes_match_on_boundaries);
836 }
837 
839  boundary_id_type this_mesh_boundary_id,
840  boundary_id_type other_mesh_boundary_id,
841  Real tol,
842  bool clear_stitched_boundary_ids,
843  bool verbose,
844  bool use_binary_search,
845  bool enforce_all_nodes_match_on_boundaries)
846 {
847  std::map<dof_id_type, dof_id_type> node_to_node_map, other_to_this_node_map; // The second is the inverse map of the first
848  std::map<dof_id_type, std::vector<dof_id_type> > node_to_elems_map;
849  // If there is only one mesh (i.e. other_mesh==NULL), then loop over this mesh twice
850  if(!other_mesh)
851  {
852  other_mesh = this;
853  }
854 
855  if( (this_mesh_boundary_id != BoundaryInfo::invalid_id) &&
856  (other_mesh_boundary_id != BoundaryInfo::invalid_id) )
857  {
858  // While finding nodes on the boundary, also find the minimum edge length
859  // of all faces on both boundaries. This will later be used in relative
860  // distance checks when stitching nodes.
862 
863  // Loop below fills in these sets for the two meshes.
864  std::set<dof_id_type> this_boundary_node_ids, other_boundary_node_ids;
865  {
866  // Make temporary fixed-size arrays for loop
867  boundary_id_type id_array[2] = {this_mesh_boundary_id, other_mesh_boundary_id};
868  std::set<dof_id_type>* set_array[2] = {&this_boundary_node_ids, &other_boundary_node_ids};
869  SerialMesh* mesh_array[2] = {this, other_mesh};
870 
871  for (unsigned i=0; i<2; ++i)
872  {
873  MeshBase::element_iterator elem_it = mesh_array[i]->elements_begin();
874  MeshBase::element_iterator elem_end = mesh_array[i]->elements_end();
875  for ( ; elem_it != elem_end; ++elem_it)
876  {
877  Elem *el = *elem_it;
878 
879  // Now check whether elem has a face on the specified boundary
880  for (unsigned int side_id=0; side_id<el->n_sides(); ++side_id)
881  if (el->neighbor(side_id) == NULL)
882  {
883  // Get *all* boundary IDs, not just the first one!
884  std::vector<boundary_id_type> bc_ids = mesh_array[i]->boundary_info->boundary_ids (el, side_id);
885 
886  if (std::count(bc_ids.begin(), bc_ids.end(), id_array[i]))
887  {
888  AutoPtr<Elem> side (el->build_side(side_id));
889  for (unsigned int node_id=0; node_id<side->n_nodes(); ++node_id)
890  set_array[i]->insert( side->node(node_id) );
891 
892  h_min = std::min(h_min, side->hmin());
893  }
894  }
895  }
896  }
897  }
898 
899  if (verbose)
900  {
901  libMesh::out << "In SerialMesh::stitch_meshes:\n"
902  << "This mesh has " << this_boundary_node_ids.size()
903  << " nodes on boundary " << this_mesh_boundary_id << ".\n"
904  << "Other mesh has " << other_boundary_node_ids.size()
905  << " nodes on boundary " << other_mesh_boundary_id << ".\n"
906  << "Minimum edge length on both surfaces is " << h_min << ".\n"
907  << std::endl;
908  }
909 
910 
911  if(use_binary_search)
912  {
913  // Store points from both stitched faces in sorted vectors for faster
914  // searching later.
915  typedef std::vector< std::pair<Point, dof_id_type> > PointVector;
916  PointVector
917  this_sorted_bndry_nodes(this_boundary_node_ids.size()),
918  other_sorted_bndry_nodes(other_boundary_node_ids.size());
919 
920  // Comparison object that will be used later. So far, I've had reasonable success
921  // with TOLERANCE...
922  FuzzyPointCompare mein_comp(TOLERANCE);
923 
924  // Create and sort the vectors we will use to do the geometric searching
925  {
926  std::set<dof_id_type>* set_array[2] = {&this_boundary_node_ids, &other_boundary_node_ids};
927  SerialMesh* mesh_array[2] = {this, other_mesh};
928  PointVector* vec_array[2] = {&this_sorted_bndry_nodes, &other_sorted_bndry_nodes};
929 
930  for (unsigned i=0; i<2; ++i)
931  {
932  std::set<dof_id_type>::iterator
933  set_it = set_array[i]->begin(),
934  set_it_end = set_array[i]->end();
935 
936  // Fill up the vector with the contents of the set...
937  for (unsigned ctr=0; set_it != set_it_end; ++set_it, ++ctr)
938  {
939  (*vec_array[i])[ctr] = std::make_pair( mesh_array[i]->point(*set_it), // The geometric point
940  *set_it ); // Its ID
941  }
942 
943  // Sort the vectors based on the FuzzyPointCompare struct op()
944  std::sort(vec_array[i]->begin(), vec_array[i]->end(), mein_comp);
945  }
946  }
947 
948  // Build up the node_to_node_map and node_to_elems_map using the sorted vectors of Points.
949  for (unsigned i=0; i<this_sorted_bndry_nodes.size(); ++i)
950  {
951  // Current point we're working on
952  Point this_point = this_sorted_bndry_nodes[i].first;
953 
954  // FuzzyPointCompare does a fuzzy equality comparison internally to handle
955  // slight differences between the list of nodes on each mesh.
956  PointVector::iterator other_iter = Utility::binary_find(other_sorted_bndry_nodes.begin(),
957  other_sorted_bndry_nodes.end(),
958  this_point,
959  mein_comp);
960 
961  // Not every node on this_sorted_bndry_nodes will necessarily be stitched, so
962  // if its pair is not found on other_mesh, just continue.
963  if (other_iter != other_sorted_bndry_nodes.end())
964  {
965  // Check that the points do indeed match - should not be necessary unless something
966  // is wrong with binary_find. To be on the safe side, we'll check.
967  {
968  // Grab the other point from the iterator
969  Point other_point = other_iter->first;
970 
971  if (!this_point.absolute_fuzzy_equals(other_point, tol*h_min))
972  {
973  libMesh::out << "Error: mismatched points: " << this_point << " and " << other_point << std::endl;
974  libmesh_error();
975  }
976  }
977 
978 
979  // Associate these two nodes in both the node_to_node_map and the other_to_this_node_map
981  this_node_id = this_sorted_bndry_nodes[i].second,
982  other_node_id = other_iter->second;
983  node_to_node_map[this_node_id] = other_node_id;
984  other_to_this_node_map[other_node_id] = this_node_id;
985  }
986 
987  }
988  }
989  else
990  {
991  // Otherwise, use a simple N^2 search to find the closest matching points. This can be helpful
992  // in the case that we have tolerance issues which cause mismatch between the two surfaces
993  // that are being stitched.
994 
995  std::set<dof_id_type>::iterator set_it = this_boundary_node_ids.begin();
996  std::set<dof_id_type>::iterator set_it_end = this_boundary_node_ids.end();
997  for( ; set_it != set_it_end; ++set_it)
998  {
999  dof_id_type this_node_id = *set_it;
1000  Node& this_node = this->node(this_node_id);
1001 
1002  bool found_matching_nodes = false;
1003 
1004  std::set<dof_id_type>::iterator other_set_it = other_boundary_node_ids.begin();
1005  std::set<dof_id_type>::iterator other_set_it_end = other_boundary_node_ids.end();
1006  for( ; other_set_it != other_set_it_end; ++other_set_it)
1007  {
1008  dof_id_type other_node_id = *other_set_it;
1009  Node& other_node = other_mesh->node(other_node_id);
1010 
1011  Real node_distance = (this_node - other_node).size();
1012 
1013  if(node_distance < tol*h_min)
1014  {
1015  // Make sure we didn't already find a matching node!
1016  if(found_matching_nodes)
1017  {
1018  libMesh::out << "Error: Found multiple matching nodes in stitch_meshes" << std::endl;
1019  libmesh_error();
1020  }
1021 
1022  node_to_node_map[this_node_id] = other_node_id;
1023  other_to_this_node_map[other_node_id] = this_node_id;
1024 
1025  found_matching_nodes = true;
1026  }
1027  }
1028  }
1029  }
1030 
1031  // Build up the node_to_elems_map, using only one loop over other_mesh
1032  {
1033  MeshBase::element_iterator other_elem_it = other_mesh->elements_begin();
1034  MeshBase::element_iterator other_elem_end = other_mesh->elements_end();
1035  for (; other_elem_it != other_elem_end; ++other_elem_it)
1036  {
1037  Elem *el = *other_elem_it;
1038 
1039  // For each node on the element, find the corresponding node
1040  // on "this" Mesh, 'this_node_id', if it exists, and push
1041  // the current element ID back onto node_to_elems_map[this_node_id].
1042  // For that we will use the reverse mapping we created at
1043  // the same time as the forward mapping.
1044  for (unsigned n=0; n<el->n_nodes(); ++n)
1045  {
1046  dof_id_type other_node_id = el->node(n);
1047  std::map<dof_id_type, dof_id_type>::iterator it =
1048  other_to_this_node_map.find(other_node_id);
1049 
1050  if (it != other_to_this_node_map.end())
1051  {
1052  dof_id_type this_node_id = it->second;
1053  node_to_elems_map[this_node_id].push_back( el->id() );
1054  }
1055  }
1056  }
1057  }
1058 
1059  if(verbose)
1060  {
1061  libMesh::out << "In SerialMesh::stitch_meshes:\n"
1062  << "Found " << node_to_node_map.size()
1063  << " matching nodes.\n"
1064  << std::endl;
1065  }
1066 
1067  if(enforce_all_nodes_match_on_boundaries)
1068  {
1069  unsigned int n_matching_nodes = node_to_node_map.size();
1070  unsigned int this_mesh_n_nodes = this_boundary_node_ids.size();
1071  unsigned int other_mesh_n_nodes = other_boundary_node_ids.size();
1072  if( (n_matching_nodes != this_mesh_n_nodes) ||
1073  (n_matching_nodes != other_mesh_n_nodes) )
1074  {
1075  libMesh::out << "Error: We expected the number of nodes to match."
1076  << std::endl;
1077  libmesh_error();
1078  }
1079  }
1080  }
1081  else
1082  {
1083  if(verbose)
1084  {
1085  libMesh::out << "Skip node merging in SerialMesh::stitch_meshes:" << std::endl;
1086  }
1087  }
1088 
1089 
1090 
1091  dof_id_type node_delta = this->n_nodes();
1092  dof_id_type elem_delta = this->n_elem();
1093 
1094  // If other_mesh!=NULL, then we have to do a bunch of work
1095  // in order to copy it to this mesh
1096  if(this!=other_mesh)
1097  {
1098  // need to increment node and element IDs of other_mesh before copying to this mesh
1099  MeshBase::node_iterator node_it = other_mesh->nodes_begin();
1100  MeshBase::node_iterator node_end = other_mesh->nodes_end();
1101  for (; node_it != node_end; ++node_it)
1102  {
1103  Node *nd = *node_it;
1104  dof_id_type new_id = nd->id() + node_delta;
1105  nd->set_id(new_id);
1106  }
1107 
1108  MeshBase::element_iterator elem_it = other_mesh->elements_begin();
1109  MeshBase::element_iterator elem_end = other_mesh->elements_end();
1110  for (; elem_it != elem_end; ++elem_it)
1111  {
1112  Elem *el = *elem_it;
1113  dof_id_type new_id = el->id() + elem_delta;
1114  el->set_id(new_id);
1115  }
1116 
1117  // Also, increment the node_to_node_map and node_to_elems_map
1118  std::map<dof_id_type, dof_id_type>::iterator node_map_it = node_to_node_map.begin();
1119  std::map<dof_id_type, dof_id_type>::iterator node_map_it_end = node_to_node_map.end();
1120  for( ; node_map_it != node_map_it_end; ++node_map_it)
1121  {
1122  node_map_it->second += node_delta;
1123  }
1124  std::map<dof_id_type, std::vector<dof_id_type> >::iterator elem_map_it = node_to_elems_map.begin();
1125  std::map<dof_id_type, std::vector<dof_id_type> >::iterator elem_map_it_end = node_to_elems_map.end();
1126  for( ; elem_map_it != elem_map_it_end; ++elem_map_it)
1127  {
1128  dof_id_type n_elems = elem_map_it->second.size();
1129  for(dof_id_type i=0; i<n_elems; i++)
1130  {
1131  (elem_map_it->second)[i] += elem_delta;
1132  }
1133  }
1134 
1135  // Copy mesh data
1136  this->copy_nodes_and_elements(*other_mesh);
1137 
1138  // Decrement node IDs of mesh to return to original state
1139  node_it = other_mesh->nodes_begin();
1140  node_end = other_mesh->nodes_end();
1141  for (; node_it != node_end; ++node_it)
1142  {
1143  Node *nd = *node_it;
1144  dof_id_type new_id = nd->id() - node_delta;
1145  nd->set_id(new_id);
1146  }
1147 
1148  elem_it = other_mesh->elements_begin();
1149  elem_end = other_mesh->elements_end();
1150  for (; elem_it != elem_end; ++elem_it)
1151  {
1152  Elem *other_elem = *elem_it;
1153 
1154  // Find the corresponding element on this mesh
1155  Elem* this_elem = this->elem(other_elem->id());
1156 
1157  // Decrement elem IDs of other_mesh to return it to original state
1158  dof_id_type new_id = other_elem->id() - elem_delta;
1159  other_elem->set_id(new_id);
1160 
1161  unsigned int n_nodes = other_elem->n_nodes();
1162  for (unsigned int n=0; n != n_nodes; ++n)
1163  {
1164  const std::vector<boundary_id_type>& ids =
1165  other_mesh->boundary_info->boundary_ids(other_elem->get_node(n));
1166  if (!ids.empty())
1167  {
1168  this->boundary_info->add_node(this_elem->get_node(n), ids);
1169  }
1170  }
1171 
1172  // Copy edge boundary info
1173  unsigned int n_edges = other_elem->n_edges();
1174  for (unsigned int edge=0; edge != n_edges; ++edge)
1175  {
1176  const std::vector<boundary_id_type>& ids =
1177  other_mesh->boundary_info->edge_boundary_ids(other_elem, edge);
1178  if (!ids.empty())
1179  {
1180  this->boundary_info->add_edge( this_elem, edge, ids);
1181  }
1182  }
1183 
1184  unsigned int n_sides = other_elem->n_sides();
1185  for (unsigned int s=0; s != n_sides; ++s)
1186  {
1187  const std::vector<boundary_id_type>& ids =
1188  other_mesh->boundary_info->boundary_ids(other_elem, s);
1189  if (!ids.empty())
1190  {
1191  this->boundary_info->add_side( this_elem, s, ids);
1192  }
1193  }
1194 
1195  }
1196 
1197  } // end if(other_mesh)
1198 
1199  // Finally, we need to "merge" the overlapping nodes
1200  // We do this by iterating over node_to_elems_map and updating
1201  // the elements so that they "point" to the nodes that came
1202  // from this mesh, rather than from other_mesh.
1203  // Then we iterate over node_to_node_map and delete the
1204  // duplicate nodes that came from other_mesh.
1205  std::map<dof_id_type, std::vector<dof_id_type> >::iterator elem_map_it = node_to_elems_map.begin();
1206  std::map<dof_id_type, std::vector<dof_id_type> >::iterator elem_map_it_end = node_to_elems_map.end();
1207  for( ; elem_map_it != elem_map_it_end; ++elem_map_it)
1208  {
1209  dof_id_type target_node_id = elem_map_it->first;
1210  dof_id_type other_node_id = node_to_node_map[target_node_id];
1211  Node& target_node = this->node(target_node_id);
1212 
1213  dof_id_type n_elems = elem_map_it->second.size();
1214  for(unsigned int i=0; i<n_elems; i++)
1215  {
1216  dof_id_type elem_id = elem_map_it->second[i];
1217  Elem* el = this->elem(elem_id);
1218 
1219  // find the local node index that we want to update
1220  unsigned int local_node_index = el->local_node(other_node_id);
1221 
1222  el->set_node(local_node_index) = &target_node;
1223  }
1224  }
1225 
1226  std::map<dof_id_type, dof_id_type>::iterator node_map_it = node_to_node_map.begin();
1227  std::map<dof_id_type, dof_id_type>::iterator node_map_it_end = node_to_node_map.end();
1228  for( ; node_map_it != node_map_it_end; ++node_map_it)
1229  {
1230  dof_id_type node_id = node_map_it->second;
1231  this->delete_node( this->node_ptr(node_id) );
1232  }
1233 
1234  this->prepare_for_use( /*skip_renumber_nodes_and_elements= */ false);
1235 
1236  // After the stitching, we may want to clear boundary IDs from element
1237  // faces that are now internal to the mesh
1238  if(clear_stitched_boundary_ids)
1239  {
1240  MeshBase::element_iterator elem_it = this->elements_begin();
1241  MeshBase::element_iterator elem_end = this->elements_end();
1242  for (; elem_it != elem_end; ++elem_it)
1243  {
1244  Elem *el = *elem_it;
1245 
1246  for (unsigned int side_id=0; side_id<el->n_sides(); side_id++)
1247  {
1248  if (el->neighbor(side_id) != NULL)
1249  {
1250  // Completely remove the side from the boundary_info object if it has either
1251  // this_mesh_boundary_id or other_mesh_boundary_id.
1252  std::vector<boundary_id_type> bc_ids = this->boundary_info->boundary_ids (el, side_id);
1253 
1254  if (std::count(bc_ids.begin(), bc_ids.end(), this_mesh_boundary_id) ||
1255  std::count(bc_ids.begin(), bc_ids.end(), other_mesh_boundary_id))
1256  this->boundary_info->remove_side(el, side_id);
1257  }
1258  }
1259  }
1260  }
1261 
1262 }
1263 
1264 
1266 {
1267  return static_cast<dof_id_type>(std::distance (this->active_elements_begin(),
1268  this->active_elements_end()));
1269 }
1270 
1271 
1272 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1274 {
1275  for (dof_id_type i=0; i<_elements.size(); ++i)
1276  if (_elements[i] && ! _elements[i]->valid_unique_id())
1277  _elements[i]->set_unique_id() = _next_unique_id++;
1278 
1279  for (dof_id_type i=0; i<_nodes.size(); ++i)
1280  if (_nodes[i] && ! _nodes[i]->valid_unique_id())
1281  _nodes[i]->set_unique_id() = _next_unique_id++;
1282 }
1283 #endif
1284 
1285 
1286 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo