mesh_tools.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 // C++ includes
21 #include <limits>
22 #include <numeric> // for std::accumulate
23 #include <set>
24 
25 // Local includes
26 #include "libmesh/elem.h"
27 #include "libmesh/elem_range.h"
28 #include "libmesh/location_maps.h"
29 #include "libmesh/mesh_base.h"
31 #include "libmesh/mesh_tools.h"
32 #include "libmesh/node_range.h"
33 #include "libmesh/parallel.h"
34 #include "libmesh/parallel_mesh.h"
35 #include "libmesh/serial_mesh.h"
36 #include "libmesh/sphere.h"
37 #include "libmesh/threads.h"
38 
39 #ifdef DEBUG
40 # include "libmesh/remote_elem.h"
41 #endif
42 
43 
44 
45 // ------------------------------------------------------------
46 // anonymous namespace for helper classes
47 namespace {
48 
49  using namespace libMesh;
50 
58  class SumElemWeight
59  {
60  public:
61  SumElemWeight () :
62  _weight(0)
63  {}
64 
65  SumElemWeight (SumElemWeight &, Threads::split) :
66  _weight(0)
67  {}
68 
69  void operator()(const ConstElemRange &range)
70  {
71  for (ConstElemRange::const_iterator it = range.begin(); it !=range.end(); ++it)
72  _weight += (*it)->n_nodes();
73  }
74 
75  dof_id_type weight() const
76  { return _weight; }
77 
78 // If we don't have threads we never need a join, and icpc yells a
79 // warning if it sees an anonymous function that's never used
80 #if LIBMESH_USING_THREADS
81  void join (const SumElemWeight &other)
82  { _weight += other.weight(); }
83 #endif
84 
85  private:
86  dof_id_type _weight;
87  };
88 
89 
96  class FindBBox
97  {
98  public:
99  FindBBox () :
100  _vmin(LIBMESH_DIM, std::numeric_limits<Real>::max()),
101  _vmax(LIBMESH_DIM, -std::numeric_limits<Real>::max())
102  {}
103 
104  FindBBox (FindBBox &other, Threads::split) :
105  _vmin(other._vmin),
106  _vmax(other._vmax)
107  {}
108 
109  std::vector<Real> & min() { return _vmin; }
110  std::vector<Real> & max() { return _vmax; }
111 
112  void operator()(const ConstNodeRange &range)
113  {
114  for (ConstNodeRange::const_iterator it = range.begin(); it != range.end(); ++it)
115  {
116  const Node *node = *it;
117  libmesh_assert(node);
118 
119  for (unsigned int i=0; i<LIBMESH_DIM; i++)
120  {
121  _vmin[i] = std::min(_vmin[i], (*node)(i));
122  _vmax[i] = std::max(_vmax[i], (*node)(i));
123  }
124  }
125  }
126 
127  void operator()(const ConstElemRange &range)
128  {
129  for (ConstElemRange::const_iterator it = range.begin(); it != range.end(); ++it)
130  {
131  const Elem *elem = *it;
132  libmesh_assert(elem);
133 
134  for (unsigned int n=0; n<elem->n_nodes(); n++)
135  {
136  const Point &point = elem->point(n);
137 
138  for (unsigned int i=0; i<LIBMESH_DIM; i++)
139  {
140  _vmin[i] = std::min(_vmin[i], point(i));
141  _vmax[i] = std::max(_vmax[i], point(i));
142  }
143  }
144  }
145  }
146 
147 // If we don't have threads we never need a join, and icpc yells a
148 // warning if it sees an anonymous function that's never used
149 #if LIBMESH_USING_THREADS
150  void join (const FindBBox &other)
151  {
152  for (unsigned int i=0; i<LIBMESH_DIM; i++)
153  {
154  _vmin[i] = std::min(_vmin[i], other._vmin[i]);
155  _vmax[i] = std::max(_vmax[i], other._vmax[i]);
156  }
157  }
158 #endif
159 
160  MeshTools::BoundingBox bbox () const
161  {
162  Point pmin(_vmin[0]
163 #if LIBMESH_DIM > 1
164  , _vmin[1]
165 #endif
166 #if LIBMESH_DIM > 2
167  , _vmin[2]
168 #endif
169  );
170  Point pmax(_vmax[0]
171 #if LIBMESH_DIM > 1
172  , _vmax[1]
173 #endif
174 #if LIBMESH_DIM > 2
175  , _vmax[2]
176 #endif
177  );
178 
179  const MeshTools::BoundingBox ret_val(pmin, pmax);
180 
181  return ret_val;
182  }
183 
184  private:
185  std::vector<Real> _vmin;
186  std::vector<Real> _vmax;
187  };
188 
189 #ifdef DEBUG
190  void assert_semiverify_dofobj(const Parallel::Communicator &communicator,
191  const DofObject *d)
192  {
193  if (d)
194  {
195  const unsigned int n_sys = d->n_systems();
196 
197  std::vector<unsigned int> n_vars (n_sys, 0);
198  for (unsigned int s = 0; s != n_sys; ++s)
199  n_vars[s] = d->n_vars(s);
200 
201  const unsigned int tot_n_vars =
202  std::accumulate(n_vars.begin(), n_vars.end(), 0);
203 
204  std::vector<unsigned int> n_comp (tot_n_vars, 0);
205  std::vector<dof_id_type> first_dof (tot_n_vars, 0);
206 
207  for (unsigned int s = 0, i=0; s != n_sys; ++s)
208  for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
209  {
210  n_comp[i] = d->n_comp(s,v);
211  first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
212  }
213 
214  libmesh_assert(communicator.semiverify(&n_sys));
215  libmesh_assert(communicator.semiverify(&n_vars));
216  libmesh_assert(communicator.semiverify(&n_comp));
217  libmesh_assert(communicator.semiverify(&first_dof));
218  }
219  else
220  {
221  const unsigned int* p_ui = NULL;
222  const std::vector<unsigned int>* p_vui = NULL;
223  const std::vector<dof_id_type>* p_vdid = NULL;
224 
225  libmesh_assert(communicator.semiverify(p_ui));
226  libmesh_assert(communicator.semiverify(p_vui));
227  libmesh_assert(communicator.semiverify(p_vui));
228  libmesh_assert(communicator.semiverify(p_vdid));
229  }
230  }
231 #endif // DEBUG
232 
233 }
234 
235 
236 namespace libMesh
237 {
238 // Small helper function to make intersect more readable.
240 {
241  return min <= check && check <= max;
242 }
243 
244 bool MeshTools::BoundingBox::intersect (const BoundingBox & other_box) const
245 {
246  // Make local variables first to make thiings more clear in a moment
247  const Real& my_min_x = this->first(0);
248  const Real& my_max_x = this->second(0);
249  const Real& other_min_x = other_box.first(0);
250  const Real& other_max_x = other_box.second(0);
251 
252  const bool x_int = is_between(my_min_x, other_min_x, my_max_x) || is_between(my_min_x, other_max_x, my_max_x) ||
253  is_between(other_min_x, my_min_x, other_max_x) || is_between(other_min_x, my_max_x, other_max_x);
254 
255  bool intersection_true = x_int;
256 
257 #if LIBMESH_DIM > 1
258  const Real& my_min_y = this->first(1);
259  const Real& my_max_y = this->second(1);
260  const Real& other_min_y = other_box.first(1);
261  const Real& other_max_y = other_box.second(1);
262 
263  const bool y_int = is_between(my_min_y, other_min_y, my_max_y) || is_between(my_min_y, other_max_y, my_max_y) ||
264  is_between(other_min_y, my_min_y, other_max_y) || is_between(other_min_y, my_max_y, other_max_y);
265 
266  intersection_true = intersection_true && y_int;
267 #endif
268 
269 #if LIBMESH_DIM > 2
270  const Real& my_min_z = this->first(2);
271  const Real& my_max_z = this->second(2);
272  const Real& other_min_z = other_box.first(2);
273  const Real& other_max_z = other_box.second(2);
274 
275  const bool z_int = is_between(my_min_z, other_min_z, my_max_z) || is_between(my_min_z, other_max_z, my_max_z) ||
276  is_between(other_min_z, my_min_z, other_max_z) || is_between(other_min_z, my_max_z, other_max_z);
277 
278  intersection_true = intersection_true && z_int;
279 #endif
280 
281  return intersection_true;
282 }
283 
285 {
286  // Make local variables first to make thiings more clear in a moment
287  Real my_min_x = this->first(0);
288  Real my_max_x = this->second(0);
289  bool x_int = is_between(my_min_x, p(0), my_max_x);
290 
291  bool intersection_true = x_int;
292 
293 #if LIBMESH_DIM > 1
294  Real my_min_y = this->first(1);
295  Real my_max_y = this->second(1);
296  bool y_int = is_between(my_min_y, p(1), my_max_y);
297 
298  intersection_true = intersection_true && y_int;
299 #endif
300 
301 
302 #if LIBMESH_DIM > 2
303  Real my_min_z = this->first(2);
304  Real my_max_z = this->second(2);
305  bool z_int = is_between(my_min_z, p(2), my_max_z);
306 
307  intersection_true = intersection_true && z_int;
308 #endif
309 
310  return intersection_true;
311 }
312 
313 // ------------------------------------------------------------
314 // MeshTools functions
316 {
317  if (!mesh.is_serial())
318  {
319  libmesh_parallel_only(mesh.comm());
321  mesh.comm().sum(weight);
322  dof_id_type unpartitioned_weight =
324  return weight + unpartitioned_weight;
325  }
326 
327  SumElemWeight sew;
328 
330  mesh.elements_end()),
331  sew);
332  return sew.weight();
333 
334 }
335 
336 
337 
339 {
340  SumElemWeight sew;
341 
343  mesh.pid_elements_end(pid)),
344  sew);
345  return sew.weight();
346 }
347 
348 
349 
351  std::vector<std::vector<dof_id_type> >& nodes_to_elem_map)
352 {
353  nodes_to_elem_map.resize (mesh.n_nodes());
354 
357 
358  for (; el != end; ++el)
359  for (unsigned int n=0; n<(*el)->n_nodes(); n++)
360  {
361  libmesh_assert_less ((*el)->node(n), nodes_to_elem_map.size());
362  libmesh_assert_less ((*el)->id(), mesh.n_elem());
363 
364  nodes_to_elem_map[(*el)->node(n)].push_back((*el)->id());
365  }
366 }
367 
368 
369 
371  std::vector<std::vector<const Elem*> >& nodes_to_elem_map)
372 {
373  nodes_to_elem_map.resize (mesh.n_nodes());
374 
377 
378  for (; el != end; ++el)
379  for (unsigned int n=0; n<(*el)->n_nodes(); n++)
380  {
381  libmesh_assert_less ((*el)->node(n), nodes_to_elem_map.size());
382 
383  nodes_to_elem_map[(*el)->node(n)].push_back(*el);
384  }
385 }
386 
387 
388 
390  std::vector<bool>& on_boundary)
391 {
392  // Resize the vector which holds boundary nodes and fill with false.
393  on_boundary.resize(mesh.n_nodes());
394  std::fill(on_boundary.begin(),
395  on_boundary.end(),
396  false);
397 
398  // Loop over elements, find those on boundary, and
399  // mark them as true in on_boundary.
402 
403  for (; el != end; ++el)
404  for (unsigned int s=0; s<(*el)->n_neighbors(); s++)
405  if ((*el)->neighbor(s) == NULL) // on the boundary
406  {
407  const AutoPtr<Elem> side((*el)->build_side(s));
408 
409  for (unsigned int n=0; n<side->n_nodes(); n++)
410  on_boundary[side->node(n)] = true;
411  }
412 }
413 
414 
415 
418 {
419  // This function must be run on all processors at once
420  libmesh_parallel_only(mesh.comm());
421 
422  FindBBox find_bbox;
423 
425  mesh.local_nodes_end()),
426  find_bbox);
427 
428  // and the unpartitioned nodes
431  find_bbox);
432 
433  // Compare the bounding boxes across processors
434  mesh.comm().min(find_bbox.min());
435  mesh.comm().max(find_bbox.max());
436 
437  return find_bbox.bbox();
438 }
439 
440 
441 
442 Sphere
444 {
445  BoundingBox bbox = bounding_box(mesh);
446 
447  const Real diag = (bbox.second - bbox.first).size();
448  const Point cent = (bbox.second + bbox.first)/2;
449 
450  return Sphere (cent, .5*diag);
451 }
452 
453 
454 
457  const processor_id_type pid)
458 {
459  libmesh_assert_less (pid, mesh.n_processors());
460 
461  FindBBox find_bbox;
462 
464  mesh.pid_elements_end(pid)),
465  find_bbox);
466 
467  return find_bbox.bbox();
468 }
469 
470 
471 
472 Sphere
474  const processor_id_type pid)
475 {
476  BoundingBox bbox = processor_bounding_box(mesh,pid);
477 
478  const Real diag = (bbox.second - bbox.first).size();
479  const Point cent = (bbox.second + bbox.first)/2;
480 
481  return Sphere (cent, .5*diag);
482 }
483 
484 
485 
488  const subdomain_id_type sid)
489 {
490  libmesh_assert_not_equal_to (mesh.n_nodes(), 0);
491 
492  Point min( 1.e30, 1.e30, 1.e30);
493  Point max(-1.e30, -1.e30, -1.e30);
494 
495  for (unsigned int e=0; e<mesh.n_elem(); e++)
496  if (mesh.elem(e)->subdomain_id() == sid)
497  for (unsigned int n=0; n<mesh.elem(e)->n_nodes(); n++)
498  for (unsigned int i=0; i<mesh.spatial_dimension(); i++)
499  {
500  min(i) = std::min(min(i), mesh.point(mesh.elem(e)->node(n))(i));
501  max(i) = std::max(max(i), mesh.point(mesh.elem(e)->node(n))(i));
502  }
503 
504  return BoundingBox (min, max);
505 }
506 
507 
508 
509 Sphere
511  const subdomain_id_type sid)
512 {
513  BoundingBox bbox = subdomain_bounding_box(mesh,sid);
514 
515  const Real diag = (bbox.second - bbox.first).size();
516  const Point cent = (bbox.second + bbox.first)/2;
517 
518  return Sphere (cent, .5*diag);
519 }
520 
521 
522 
524  std::vector<ElemType>& et)
525 {
528 
529  // Automatically get the first type
530  et.push_back((*el)->type()); ++el;
531 
532  // Loop over the rest of the elements.
533  // If the current element type isn't in the
534  // vector, insert it.
535  for (; el != end; ++el)
536  if (!std::count(et.begin(), et.end(), (*el)->type()))
537  et.push_back((*el)->type());
538 }
539 
540 
541 
543  const ElemType type)
544 {
545  return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
546  mesh.type_elements_end (type)));
547 }
548 
549 
550 
552  const ElemType type)
553 {
554  return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
555  mesh.active_type_elements_end (type)));
556 }
557 
559  const ElemType type,
560  const unsigned int level)
561 {
562  dof_id_type cnt = 0;
563  // iterate over the elements of the specified type
566 
567  for(; el!=end; ++el)
568  if( ((*el)->level() == level) && !(*el)->subactive())
569  cnt++;
570 
571  return cnt;
572 }
573 
574 
576 {
577  unsigned int max_level = 0;
578 
581 
582  for( ; el != end_el; ++el)
583  max_level = std::max((*el)->level(), max_level);
584 
585  return max_level + 1;
586 }
587 
588 
589 
591 {
592  libmesh_parallel_only(mesh.comm());
593 
594  unsigned int nl = MeshTools::n_active_local_levels(mesh);
595 
598  const MeshBase::const_element_iterator end_el =
600 
601  for( ; el != end_el; ++el)
602  if ((*el)->active())
603  nl = std::max((*el)->level() + 1, nl);
604 
605  mesh.comm().max(nl);
606  return nl;
607 }
608 
609 
610 
612 {
613  unsigned int max_level = 0;
614 
617 
618  for( ; el != end_el; ++el)
619  max_level = std::max((*el)->level(), max_level);
620 
621  return max_level + 1;
622 }
623 
624 
625 
626 unsigned int MeshTools::n_levels(const MeshBase& mesh)
627 {
628  libmesh_parallel_only(mesh.comm());
629 
630  unsigned int nl = MeshTools::n_local_levels(mesh);
631 
634  const MeshBase::const_element_iterator end_el =
636 
637  for( ; el != end_el; ++el)
638  nl = std::max((*el)->level() + 1, nl);
639 
640  mesh.comm().max(nl);
641  return nl;
642 }
643 
644 
645 
647  std::set<dof_id_type>& not_subactive_node_ids)
648 {
650  const MeshBase::const_element_iterator end_el = mesh.elements_end();
651  for( ; el != end_el; ++el)
652  {
653  const Elem* elem = (*el);
654  if(!elem->subactive())
655  for (unsigned int n=0; n<elem->n_nodes(); ++n)
656  not_subactive_node_ids.insert(elem->node(n));
657  }
658 }
659 
660 
661 
664 {
665  return std::distance(begin, end);
666 }
667 
668 
669 
672 {
673  return std::distance(begin, end);
674 }
675 
676 
677 
678 unsigned int MeshTools::n_p_levels (const MeshBase& mesh)
679 {
680  libmesh_parallel_only(mesh.comm());
681 
682  unsigned int max_p_level = 0;
683 
684  // first my local elements
686  el = mesh.local_elements_begin(),
687  end_el = mesh.local_elements_end();
688 
689  for( ; el != end_el; ++el)
690  max_p_level = std::max((*el)->p_level(), max_p_level);
691 
692  // then any unpartitioned objects
693  el = mesh.unpartitioned_elements_begin();
694  end_el = mesh.unpartitioned_elements_end();
695 
696  for( ; el != end_el; ++el)
697  max_p_level = std::max((*el)->p_level(), max_p_level);
698 
699  mesh.comm().max(max_p_level);
700  return max_p_level + 1;
701 }
702 
703 
704 
706  std::vector<std::vector<const Elem*> >& nodes_to_elem_map,
707  std::vector<const Node*>& neighbors)
708 {
709  dof_id_type global_id = n.id();
710 
711  //Iterators to iterate through the elements that include this node
712  std::vector<const Elem*>::const_iterator el = nodes_to_elem_map[global_id].begin();
713  std::vector<const Elem*>::const_iterator end_el = nodes_to_elem_map[global_id].end();
714 
715  unsigned int n_ed=0; // Number of edges on the element
716  unsigned int ed=0; // Current edge
717  unsigned int l_n=0; // Local node number
718  unsigned int o_n=0; // Other node on this edge
719 
720  //Assume we find a edge... then prove ourselves wrong...
721  bool found_edge=true;
722 
723  Node * node_to_save = NULL;
724 
725  //Look through the elements that contain this node
726  //find the local node id... then find the side that
727  //node lives on in the element
728  //next, look for the _other_ node on that side
729  //That other node is a "nodal_neighbor"... save it
730  for(;el != end_el;el++)
731  {
732  //We only care about active elements...
733  if((*el)->active())
734  {
735  n_ed=(*el)->n_edges();
736 
737  //Find the local node id
738  while(global_id != (*el)->node(l_n++)) { }
739  l_n--; //Hmmm... take the last one back off
740 
741  while(ed<n_ed)
742  {
743 
744  //Find the edge the node is on
745  while(found_edge && !(*el)->is_node_on_edge(l_n,ed++))
746  {
747  //This only happens if all the edges have already been found
748  if(ed>=n_ed)
749  found_edge=false;
750  }
751 
752  //Did we find one?
753  if(found_edge)
754  {
755  ed--; //Take the last one back off again
756 
757  //Now find the other node on that edge
758  while(!(*el)->is_node_on_edge(o_n++,ed) || global_id==(*el)->node(o_n-1)) { }
759  o_n--;
760 
761  //We've found one! Save it..
762  node_to_save=(*el)->get_node(o_n);
763 
764  //Search to see if we've already found this one
765  std::vector<const Node*>::const_iterator result = std::find(neighbors.begin(),neighbors.end(),node_to_save);
766 
767  //If we didn't find it and add it to the vector
768  if(result == neighbors.end())
769  neighbors.push_back(node_to_save);
770  }
771 
772  //Reset to look for another
773  o_n=0;
774 
775  //Keep looking for edges, node may be on more than one edge
776  ed++;
777  }
778 
779  //Reset to get ready for the next element
780  l_n=ed=0;
781  found_edge=true;
782  }
783  }
784 }
785 
786 void MeshTools::find_hanging_nodes_and_parents(const MeshBase& mesh, std::map<dof_id_type, std::vector<dof_id_type> >& hanging_nodes)
787 {
790 
791  //Loop through all the elements
792  for (; it != end; ++it)
793  {
794  //Save it off for easier access
795  const Elem* elem = (*it);
796 
797  //Right now this only works for quad4's
798  //libmesh_assert_equal_to (elem->type(), libMeshEnums::QUAD4);
799  if(elem->type() == libMeshEnums::QUAD4)
800  {
801  //Loop over the sides looking for sides that have hanging nodes
802  //This code is inspired by compute_proj_constraints()
803  for (unsigned int s=0; s<elem->n_sides(); s++)
804  {
805  //If not a boundary node
806  if (elem->neighbor(s) != NULL)
807  {
808  // Get pointers to the element's neighbor.
809  const Elem* neigh = elem->neighbor(s);
810 
811  //Is there a coarser element next to this one?
812  if (neigh->level() < elem->level())
813  {
814  const Elem *ancestor = elem;
815  while (neigh->level() < ancestor->level())
816  ancestor = ancestor->parent();
817  unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
818  libmesh_assert_less (s_neigh, neigh->n_neighbors());
819 
820  //Couple of helper uints...
821  unsigned int local_node1=0;
822  unsigned int local_node2=0;
823 
824  bool found_in_neighbor = false;
825 
826  //Find the two vertices that make up this side
827  while(!elem->is_node_on_side(local_node1++,s)) { }
828  local_node1--;
829 
830  //Start looking for the second one with the next node
831  local_node2=local_node1+1;
832 
833  //Find the other one
834  while(!elem->is_node_on_side(local_node2++,s)) { }
835  local_node2--;
836 
837  //Pull out their global ids:
838  dof_id_type node1 = elem->node(local_node1);
839  dof_id_type node2 = elem->node(local_node2);
840 
841  //Now find which node is present in the neighbor
842  //FIXME This assumes a level one rule!
843  //The _other_ one is the hanging node
844 
845  //First look for the first one
846  //FIXME could be streamlined a bit
847  for(unsigned int n=0;n<neigh->n_sides();n++)
848  {
849  if(neigh->node(n) == node1)
850  found_in_neighbor=true;
851  }
852 
853  dof_id_type hanging_node=0;
854 
855  if(!found_in_neighbor)
856  hanging_node=node1;
857  else //If it wasn't node1 then it must be node2!
858  hanging_node=node2;
859 
860  //Reset these for reuse
861  local_node1=0;
862  local_node2=0;
863 
864  //Find the first node that makes up the side in the neighbor (these should be the parent nodes)
865  while(!neigh->is_node_on_side(local_node1++,s_neigh)) { }
866  local_node1--;
867 
868  local_node2=local_node1+1;
869 
870  //Find the second node...
871  while(!neigh->is_node_on_side(local_node2++,s_neigh)) { }
872  local_node2--;
873 
874  //Save them if we haven't already found the parents for this one
875  if(hanging_nodes[hanging_node].size()<2)
876  {
877  hanging_nodes[hanging_node].push_back(neigh->node(local_node1));
878  hanging_nodes[hanging_node].push_back(neigh->node(local_node2));
879  }
880  }
881  }
882  }
883  }
884  }
885 }
886 
887 
888 
891  LocationMap<Node> &loc_map)
892 {
893  // This function must be run on all processors at once
894  libmesh_parallel_only(mesh.comm());
895 
896  // We'll need the new_nodes_map to answer other processors'
897  // requests. It should never be empty unless we don't have any
898  // nodes.
899  libmesh_assert(mesh.nodes_begin() == mesh.nodes_end() ||
900  !loc_map.empty());
901 
902  // Fix all nodes' processor ids. Coarsening may have left us with
903  // nodes which are no longer touched by any elements of the same
904  // processor id, and for DofMap to work we need to fix that.
905 
906  // In the first pass, invalidate processor ids for nodes on active
907  // elements. We avoid touching subactive-only nodes.
909  const MeshBase::element_iterator e_end = mesh.active_elements_end();
910  for (; e_it != e_end; ++e_it)
911  {
912  Elem *elem = *e_it;
913  for (unsigned int n=0; n != elem->n_nodes(); ++n)
914  {
915  Node *node = elem->get_node(n);
916  node->invalidate_processor_id();
917  }
918  }
919 
920  // In the second pass, find the lowest processor ids on active
921  // elements touching each node, and set the node processor id.
922  for (e_it = mesh.active_elements_begin(); e_it != e_end; ++e_it)
923  {
924  Elem *elem = *e_it;
925  processor_id_type proc_id = elem->processor_id();
926  for (unsigned int n=0; n != elem->n_nodes(); ++n)
927  {
928  Node *node = elem->get_node(n);
930  node->processor_id() > proc_id)
931  node->processor_id() = proc_id;
932  }
933  }
934 
935  // Those two passes will correct every node that touches a local
936  // element, but we can't be sure about nodes touching remote
937  // elements. Fix those now.
939  (mesh, loc_map);
940 }
941 
942 
943 
944 #ifdef DEBUG
946 {
948  mesh.elements_begin();
949  const MeshBase::const_element_iterator el_end =
950  mesh.elements_end();
951  if (el == el_end)
952  return;
953 
954  const unsigned int n_sys = (*el)->n_systems();
955 
956  for (; el != el_end; ++el)
957  {
958  const Elem *elem = *el;
959  libmesh_assert_equal_to (elem->n_systems(), n_sys);
960  }
961 
963  mesh.nodes_begin();
964  const MeshBase::const_node_iterator node_end =
965  mesh.nodes_end();
966 
967  if (node_it == node_end)
968  return;
969 
970  for (; node_it != node_end; ++node_it)
971  {
972  const Node *node = *node_it;
973  libmesh_assert_equal_to (node->n_systems(), n_sys);
974  }
975 }
976 
977 
978 
980 {
981 #ifdef LIBMESH_ENABLE_AMR
983  mesh.elements_begin();
984  const MeshBase::const_element_iterator el_end =
985  mesh.elements_end();
986 
987  for (; el != el_end; ++el)
988  {
989  const Elem *elem = *el;
990 
991  if (elem->refinement_flag() == Elem::JUST_REFINED ||
992  elem->refinement_flag() == Elem::INACTIVE)
993  continue;
994 
995  if (elem->has_dofs())
997 
998  for (unsigned int n=0; n != elem->n_nodes(); ++n)
999  {
1000  const Node *node = elem->get_node(n);
1001  if (node->has_dofs())
1003  }
1004  }
1005 #endif // LIBMESH_ENABLE_AMR
1006 }
1007 
1008 
1009 
1011 {
1012  const MeshBase::const_element_iterator el_end =
1013  mesh.elements_end();
1015  mesh.elements_begin(); el != el_end; ++el)
1016  {
1017  const Elem* elem = *el;
1018  libmesh_assert (elem);
1019  while (elem)
1020  {
1022  for (unsigned int n=0; n != elem->n_neighbors(); ++n)
1023  if (elem->neighbor(n) &&
1024  elem->neighbor(n) != remote_elem)
1026 
1027  libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1028  elem = elem->parent();
1029  }
1030  }
1031 }
1032 
1033 
1035 {
1036  const MeshBase::const_element_iterator el_end =
1037  mesh.local_elements_end();
1039  mesh.local_elements_begin(); el != el_end; ++el)
1040  {
1041  const Elem* elem = *el;
1042  libmesh_assert (elem);
1043  for (unsigned int n=0; n != elem->n_neighbors(); ++n)
1044  libmesh_assert_not_equal_to (elem->neighbor(n), remote_elem);
1045 #ifdef LIBMESH_ENABLE_AMR
1046  const Elem* parent = elem->parent();
1047  if (parent)
1048  {
1049  libmesh_assert_not_equal_to (parent, remote_elem);
1050  for (unsigned int c=0; c != elem->n_children(); ++c)
1051  libmesh_assert_not_equal_to (parent->child(c), remote_elem);
1052  }
1053 #endif
1054  }
1055 }
1056 
1057 
1059  const Elem *bad_elem)
1060 {
1061  const MeshBase::const_element_iterator el_end =
1062  mesh.elements_end();
1064  mesh.elements_begin(); el != el_end; ++el)
1065  {
1066  const Elem* elem = *el;
1067  libmesh_assert (elem);
1068  libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1069  for (unsigned int n=0; n != elem->n_neighbors(); ++n)
1070  libmesh_assert_not_equal_to (elem->neighbor(n), bad_elem);
1071 #ifdef LIBMESH_ENABLE_AMR
1072  if (elem->has_children())
1073  for (unsigned int c=0; c != elem->n_children(); ++c)
1074  libmesh_assert_not_equal_to (elem->child(c), bad_elem);
1075 #endif
1076  }
1077 }
1078 
1079 
1080 
1082 {
1083  processor_id_type lastprocid = 0;
1084  dof_id_type lastelemid = 0;
1085 
1086  const MeshBase::const_element_iterator el_end =
1087  mesh.active_elements_end();
1089  mesh.active_elements_begin(); el != el_end; ++el)
1090  {
1091  const Elem* elem = *el;
1092  libmesh_assert (elem);
1093  processor_id_type elemprocid = elem->processor_id();
1094  dof_id_type elemid = elem->id();
1095 
1096  libmesh_assert_greater_equal (elemid, lastelemid);
1097  libmesh_assert_greater_equal (elemprocid, lastprocid);
1098 
1099  lastelemid = elemid;
1100  lastprocid = elemprocid;
1101  }
1102 }
1103 
1104 
1105 
1107 {
1108  const MeshBase::const_element_iterator el_end =
1109  mesh.elements_end();
1111  mesh.elements_begin(); el != el_end; ++el)
1112  {
1113  const Elem* elem = *el;
1114  libmesh_assert (elem);
1115 
1116  const Elem* parent = elem->parent();
1117 
1118  if (parent)
1119  {
1120  libmesh_assert_greater_equal (elem->id(), parent->id());
1121  libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1122  }
1123  }
1124 }
1125 
1126 
1127 
1129 {
1130  std::set<const Node*> used_nodes;
1131 
1132  const MeshBase::const_element_iterator el_end =
1133  mesh.elements_end();
1135  mesh.elements_begin(); el != el_end; ++el)
1136  {
1137  const Elem* elem = *el;
1138  libmesh_assert (elem);
1139 
1140  for (unsigned int n=0; n<elem->n_nodes(); ++n)
1141  used_nodes.insert(elem->get_node(n));
1142  }
1143 
1144  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
1145 
1146  for (MeshBase::const_node_iterator node_it = mesh.nodes_begin();
1147  node_it != node_end; ++node_it)
1148  {
1149  Node *node = *node_it;
1150  libmesh_assert(node);
1151  libmesh_assert(used_nodes.count(node));
1152  }
1153 }
1154 
1155 
1156 
1157 namespace MeshTools {
1158 
1160 {
1161  if (mesh.n_processors() == 1)
1162  return;
1163 
1164  libmesh_parallel_only(mesh.comm());
1165 
1166  dof_id_type pmax_elem_id = mesh.max_elem_id();
1167  mesh.comm().max(pmax_elem_id);
1168 
1169  for (dof_id_type i=0; i != pmax_elem_id; ++i)
1170  assert_semiverify_dofobj(mesh.comm(),
1171  mesh.query_elem(i));
1172 
1173  dof_id_type pmax_node_id = mesh.max_node_id();
1174  mesh.comm().max(pmax_node_id);
1175 
1176  for (dof_id_type i=0; i != pmax_node_id; ++i)
1177  assert_semiverify_dofobj(mesh.comm(),
1178  mesh.query_node_ptr(i));
1179 }
1180 
1181 template <>
1183 {
1184  if (mesh.n_processors() == 1)
1185  return;
1186 
1187  libmesh_parallel_only(mesh.comm());
1188 
1189  // We want this test to be valid even when called even after nodes
1190  // have been added asynchonously but before they're renumbered
1191  dof_id_type parallel_max_elem_id = mesh.max_elem_id();
1192  mesh.comm().max(parallel_max_elem_id);
1193 
1194  // Check processor ids for consistency between processors
1195 
1196  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
1197  {
1198  const Elem *elem = mesh.query_elem(i);
1199 
1200  processor_id_type min_id =
1201  elem ? elem->processor_id() :
1203  mesh.comm().min(min_id);
1204 
1205  processor_id_type max_id =
1206  elem ? elem->processor_id() :
1208  mesh.comm().max(max_id);
1209 
1210  if (elem)
1211  {
1212  libmesh_assert_equal_to (min_id, elem->processor_id());
1213  libmesh_assert_equal_to (max_id, elem->processor_id());
1214  }
1215 
1216  if (min_id == mesh.processor_id())
1217  libmesh_assert(elem);
1218  }
1219 
1220  // If we're adaptively refining, check processor ids for consistency
1221  // between parents and children.
1222 #ifdef LIBMESH_ENABLE_AMR
1223 
1224  // Ancestor elements we won't worry about, but subactive and active
1225  // elements ought to have parents with consistent processor ids
1226 
1227  const MeshBase::const_element_iterator el_end =
1228  mesh.elements_end();
1230  mesh.elements_begin(); el != el_end; ++el)
1231  {
1232  const Elem *elem = *el;
1233  libmesh_assert(elem);
1234 
1235  if (!elem->active() && !elem->subactive())
1236  continue;
1237 
1238  const Elem *parent = elem->parent();
1239 
1240  if (parent)
1241  {
1242  libmesh_assert(parent->has_children());
1243  processor_id_type parent_procid = parent->processor_id();
1244  bool matching_child_id = false;
1245  for (unsigned int c = 0; c != parent->n_children(); ++c)
1246  {
1247  const Elem* child = parent->child(c);
1248  libmesh_assert(child);
1249 
1250  // If we've got a remote_elem then we don't know whether
1251  // it's responsible for the parent's processor id; all
1252  // we can do is assume it is and let its processor fail
1253  // an assert if there's something wrong.
1254  if (child == remote_elem ||
1255  child->processor_id() == parent_procid)
1256  matching_child_id = true;
1257  }
1258  libmesh_assert(matching_child_id);
1259  }
1260  }
1261 #endif
1262 }
1263 
1264 
1265 
1266 template <>
1268 {
1269  if (mesh.n_processors() == 1)
1270  return;
1271 
1272  libmesh_parallel_only(mesh.comm());
1273 
1274  // We want this test to be valid even when called even after nodes
1275  // have been added asynchonously but before they're renumbered
1276  dof_id_type parallel_max_node_id = mesh.max_node_id();
1277  mesh.comm().max(parallel_max_node_id);
1278 
1279  // Check processor ids for consistency between processors
1280 
1281  for (dof_id_type i=0; i != parallel_max_node_id; ++i)
1282  {
1283  const Node *node = mesh.query_node_ptr(i);
1284 
1285  processor_id_type min_id =
1286  node ? node->processor_id() :
1288  mesh.comm().min(min_id);
1289 
1290  processor_id_type max_id =
1291  node ? node->processor_id() :
1293  mesh.comm().max(max_id);
1294 
1295  if (node)
1296  {
1297  libmesh_assert_equal_to (min_id, node->processor_id());
1298  libmesh_assert_equal_to (max_id, node->processor_id());
1299  }
1300 
1301  if (min_id == mesh.processor_id())
1302  libmesh_assert(node);
1303  }
1304 
1305  std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1306 
1307  const MeshBase::const_element_iterator el_end =
1310  mesh.active_local_elements_begin(); el != el_end; ++el)
1311  {
1312  const Elem* elem = *el;
1313  libmesh_assert (elem);
1314 
1315  for (unsigned int i=0; i != elem->n_nodes(); ++i)
1316  {
1317  const Node *node = elem->get_node(i);
1318  dof_id_type nodeid = node->id();
1319  node_touched_by_me[nodeid] = true;
1320  }
1321  }
1322  std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1323  mesh.comm().max(node_touched_by_anyone);
1324 
1325  const MeshBase::const_node_iterator nd_end = mesh.local_nodes_end();
1327  nd != nd_end; ++nd)
1328  {
1329  const Node *node = *nd;
1330  libmesh_assert(node);
1331 
1332  dof_id_type nodeid = node->id();
1333  libmesh_assert(!node_touched_by_anyone[nodeid] ||
1334  node_touched_by_me[nodeid]);
1335  }
1336 }
1337 
1338 } // namespace MeshTools
1339 
1340 
1341 
1342 #ifdef LIBMESH_ENABLE_AMR
1344 {
1345  libmesh_parallel_only(mesh.comm());
1346  if (mesh.n_processors() == 1)
1347  return;
1348 
1349  std::vector<unsigned char> my_elem_h_state(mesh.max_elem_id(), 255);
1350  std::vector<unsigned char> my_elem_p_state(mesh.max_elem_id(), 255);
1351 
1352  const MeshBase::const_element_iterator el_end =
1353  mesh.elements_end();
1355  mesh.elements_begin(); el != el_end; ++el)
1356  {
1357  const Elem* elem = *el;
1358  libmesh_assert (elem);
1359  dof_id_type elemid = elem->id();
1360 
1361  my_elem_h_state[elemid] =
1362  static_cast<unsigned char>(elem->refinement_flag());
1363 
1364  my_elem_p_state[elemid] =
1365  static_cast<unsigned char>(elem->p_refinement_flag());
1366  }
1367  std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
1368  mesh.comm().min(min_elem_h_state);
1369 
1370  std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
1371  mesh.comm().min(min_elem_p_state);
1372 
1373  for (dof_id_type i=0; i!= mesh.max_elem_id(); ++i)
1374  {
1375  libmesh_assert(my_elem_h_state[i] == 255 ||
1376  my_elem_h_state[i] == min_elem_h_state[i]);
1377  libmesh_assert(my_elem_p_state[i] == 255 ||
1378  my_elem_p_state[i] == min_elem_p_state[i]);
1379  }
1380 }
1381 #else
1383 {
1384 }
1385 #endif // LIBMESH_ENABLE_AMR
1386 
1387 
1388 
1389 #ifdef LIBMESH_ENABLE_AMR
1391 {
1392  const MeshBase::const_element_iterator el_end =
1393  mesh.elements_end();
1395  mesh.elements_begin(); el != el_end; ++el)
1396  {
1397  const Elem *elem = *el;
1398  libmesh_assert(elem);
1399  if (elem->has_children())
1400  for (unsigned int n=0; n != elem->n_children(); ++n)
1401  {
1402  libmesh_assert(elem->child(n));
1403  if (elem->child(n) != remote_elem)
1404  libmesh_assert_equal_to (elem->child(n)->parent(), elem);
1405  }
1406  if (elem->active())
1407  {
1408  libmesh_assert(!elem->ancestor());
1409  libmesh_assert(!elem->subactive());
1410  }
1411  else if (elem->ancestor())
1412  {
1413  libmesh_assert(!elem->subactive());
1414  }
1415  else
1416  libmesh_assert(elem->subactive());
1417  }
1418 }
1419 #else
1421 {
1422 }
1423 #endif // LIBMESH_ENABLE_AMR
1424 
1425 
1426 
1428 {
1429  const MeshBase::const_element_iterator el_end = mesh.elements_end();
1431  el != el_end; ++el)
1432  {
1433  const Elem* elem = *el;
1434  libmesh_assert (elem);
1436  }
1437 }
1438 
1439 
1440 
1441 #endif
1442 
1443 
1444 
1446 {
1448 }
1449 
1450 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo