elem.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 <algorithm> // for std::sort
22 #include <iterator> // for std::ostream_iterator
23 #include <sstream>
24 #include <limits> // for std::numeric_limits<>
25 #include <cmath> // for std::sqrt()
26 
27 // Local includes
28 #include "libmesh/elem.h"
29 #include "libmesh/fe_type.h"
30 #include "libmesh/fe_interface.h"
31 #include "libmesh/edge_edge2.h"
32 #include "libmesh/edge_edge3.h"
33 #include "libmesh/edge_edge4.h"
34 #include "libmesh/edge_inf_edge2.h"
35 #include "libmesh/face_tri3.h"
36 #include "libmesh/face_tri6.h"
37 #include "libmesh/face_quad4.h"
38 #include "libmesh/face_quad8.h"
39 #include "libmesh/face_quad9.h"
40 #include "libmesh/face_inf_quad4.h"
41 #include "libmesh/face_inf_quad6.h"
42 #include "libmesh/cell_tet4.h"
43 #include "libmesh/cell_tet10.h"
44 #include "libmesh/cell_hex8.h"
45 #include "libmesh/cell_hex20.h"
46 #include "libmesh/cell_hex27.h"
47 #include "libmesh/cell_inf_hex8.h"
48 #include "libmesh/cell_inf_hex16.h"
49 #include "libmesh/cell_inf_hex18.h"
50 #include "libmesh/cell_prism6.h"
51 #include "libmesh/cell_prism15.h"
52 #include "libmesh/cell_prism18.h"
55 #include "libmesh/cell_pyramid5.h"
56 #include "libmesh/cell_pyramid14.h"
57 #include "libmesh/fe_base.h"
58 #include "libmesh/mesh_base.h"
60 #include "libmesh/remote_elem.h"
61 #include "libmesh/reference_elem.h"
62 #include "libmesh/string_to_enum.h"
63 
64 #ifdef LIBMESH_ENABLE_PERIODIC
65 #include "libmesh/mesh.h"
67 #include "libmesh/boundary_info.h"
68 #endif
69 
70 namespace libMesh
71 {
72 
73 // Initialize static member variables
74 const unsigned int Elem::type_to_n_nodes_map [] =
75  {
76  2, // EDGE2
77  3, // EDGE3
78  4, // EDGE4
79 
80  3, // TRI3
81  6, // TRI6
82 
83  4, // QUAD4
84  8, // QUAD8
85  9, // QUAD9
86 
87  4, // TET4
88  10, // TET10
89 
90  8, // HEX8
91  20, // HEX20
92  27, // HEX27
93 
94  6, // PRISM6
95  15, // PRISM15
96  18, // PRISM18
97 
98  5, // PYRAMID5
99  14, // PYRAMID14
100 
101  2, // INFEDGE2
102 
103  4, // INFQUAD4
104  6, // INFQUAD6
105 
106  8, // INFHEX8
107  16, // INFHEX16
108  18, // INFHEX18
109 
110  6, // INFPRISM6
111  16, // INFPRISM12
112 
113  1, // NODEELEM
114  };
115 
116 const unsigned int Elem::type_to_n_sides_map [] =
117  {
118  2, // EDGE2
119  2, // EDGE3
120  2, // EDGE4
121 
122  3, // TRI3
123  3, // TRI6
124 
125  4, // QUAD4
126  4, // QUAD8
127  4, // QUAD9
128 
129  4, // TET4
130  4, // TET10
131 
132  6, // HEX8
133  6, // HEX20
134  6, // HEX27
135 
136  5, // PRISM6
137  5, // PRISM15
138  5, // PRISM18
139 
140  5, // PYRAMID5
141  5, // PYRAMID14
142 
143  2, // INFEDGE2
144 
145  3, // INFQUAD4
146  3, // INFQUAD6
147 
148  5, // INFHEX8
149  5, // INFHEX16
150  5, // INFHEX18
151 
152  4, // INFPRISM6
153  4, // INFPRISM12
154 
155  0, // NODEELEM
156  };
157 
158 const unsigned int Elem::type_to_n_edges_map [] =
159  {
160  0, // EDGE2
161  0, // EDGE3
162  0, // EDGE4
163 
164  3, // TRI3
165  3, // TRI6
166 
167  4, // QUAD4
168  4, // QUAD8
169  4, // QUAD9
170 
171  6, // TET4
172  6, // TET10
173 
174  12, // HEX8
175  12, // HEX20
176  12, // HEX27
177 
178  9, // PRISM6
179  9, // PRISM15
180  9, // PRISM18
181 
182  8, // PYRAMID5
183  8, // PYRAMID14
184 
185  0, // INFEDGE2
186 
187  4, // INFQUAD4
188  4, // INFQUAD6
189 
190  8, // INFHEX8
191  8, // INFHEX16
192  8, // INFHEX18
193 
194  6, // INFPRISM6
195  6, // INFPRISM12
196 
197  0, // NODEELEM
198  };
199 
200 // ------------------------------------------------------------
201 // Elem class member funcions
203  Elem* p)
204 {
205  Elem* elem = NULL;
206 
207  switch (type)
208  {
209  // 1D elements
210  case EDGE2:
211  {
212  elem = new Edge2(p);
213  break;
214  }
215  case EDGE3:
216  {
217  elem = new Edge3(p);
218  break;
219  }
220  case EDGE4:
221  {
222  elem = new Edge4(p);
223  break;
224  }
225 
226 
227 
228  // 2D elements
229  case TRI3:
230  {
231  elem = new Tri3(p);
232  break;
233  }
234  case TRI6:
235  {
236  elem = new Tri6(p);
237  break;
238  }
239  case QUAD4:
240  {
241  elem = new Quad4(p);
242  break;
243  }
244  case QUAD8:
245  {
246  elem = new Quad8(p);
247  break;
248  }
249  case QUAD9:
250  {
251  elem = new Quad9(p);
252  break;
253  }
254 
255 
256  // 3D elements
257  case TET4:
258  {
259  elem = new Tet4(p);
260  break;
261  }
262  case TET10:
263  {
264  elem = new Tet10(p);
265  break;
266  }
267  case HEX8:
268  {
269  elem = new Hex8(p);
270  break;
271  }
272  case HEX20:
273  {
274  elem = new Hex20(p);
275  break;
276  }
277  case HEX27:
278  {
279  elem = new Hex27(p);
280  break;
281  }
282  case PRISM6:
283  {
284  elem = new Prism6(p);
285  break;
286  }
287  case PRISM15:
288  {
289  elem = new Prism15(p);
290  break;
291  }
292  case PRISM18:
293  {
294  elem = new Prism18(p);
295  break;
296  }
297  case PYRAMID5:
298  {
299  elem = new Pyramid5(p);
300  break;
301  }
302  case PYRAMID14:
303  {
304  elem = new Pyramid14(p);
305  break;
306  }
307 
308 
309 
310 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
311 
312  // 1D infinite elements
313  case INFEDGE2:
314  {
315  elem = new InfEdge2(p);
316  break;
317  }
318 
319 
320  // 2D infinite elements
321  case INFQUAD4:
322  {
323  elem = new InfQuad4(p);
324  break;
325  }
326  case INFQUAD6:
327  {
328  elem = new InfQuad6(p);
329  break;
330  }
331 
332 
333  // 3D infinite elements
334  case INFHEX8:
335  {
336  elem = new InfHex8(p);
337  break;
338  }
339  case INFHEX16:
340  {
341  elem = new InfHex16(p);
342  break;
343  }
344  case INFHEX18:
345  {
346  elem = new InfHex18(p);
347  break;
348  }
349  case INFPRISM6:
350  {
351  elem = new InfPrism6(p);
352  break;
353  }
354  case INFPRISM12:
355  {
356  elem = new InfPrism12(p);
357  break;
358  }
359 
360 #endif
361 
362  default:
363  {
364  libMesh::err << "ERROR: Undefined element type!." << std::endl;
365  libmesh_error();
366  }
367  }
368 
369 
370  AutoPtr<Elem> ap(elem);
371  return ap;
372 }
373 
374 
375 
376 const Elem* Elem::reference_elem () const
377 {
378  return &(ReferenceElem::get(this->type()));
379 }
380 
381 
382 
384 {
385  Point cp;
386 
387  for (unsigned int n=0; n<this->n_vertices(); n++)
388  cp.add (this->point(n));
389 
390  return (cp /= static_cast<Real>(this->n_vertices()));
391 }
392 
393 
394 
396 {
398 
399  for (unsigned int n_outer=0; n_outer<this->n_vertices(); n_outer++)
400  for (unsigned int n_inner=n_outer+1; n_inner<this->n_vertices(); n_inner++)
401  {
402  const Point diff = (this->point(n_outer) - this->point(n_inner));
403 
404  h_min = std::min(h_min,diff.size_sq());
405  }
406 
407  return std::sqrt(h_min);
408 }
409 
410 
411 
413 {
414  Real h_max=0;
415 
416  for (unsigned int n_outer=0; n_outer<this->n_vertices(); n_outer++)
417  for (unsigned int n_inner=n_outer+1; n_inner<this->n_vertices(); n_inner++)
418  {
419  const Point diff = (this->point(n_outer) - this->point(n_inner));
420 
421  h_max = std::max(h_max,diff.size_sq());
422  }
423 
424  return std::sqrt(h_max);
425 }
426 
427 
428 
429 Real Elem::length(const unsigned int n1,
430  const unsigned int n2) const
431 {
432  libmesh_assert_less ( n1, this->n_vertices() );
433  libmesh_assert_less ( n2, this->n_vertices() );
434 
435  return (this->point(n1) - this->point(n2)).size();
436 }
437 
438 
439 
440 bool Elem::operator == (const Elem& rhs) const
441 {
442 
443  // Cast rhs to an Elem*
444 // const Elem* rhs_elem = dynamic_cast<const Elem*>(&rhs);
445  const Elem* rhs_elem = &rhs;
446 
447  // If we cannot cast to an Elem*, rhs must be a Node
448 // if(rhs_elem == static_cast<const Elem*>(NULL))
449 // return false;
450 
451 // libmesh_assert (n_nodes());
452 // libmesh_assert (rhs.n_nodes());
453 
454 // // Elements can only be equal if they
455 // // contain the same number of nodes.
456 // if (this->n_nodes() == rhs.n_nodes())
457 // {
458 // // Create a set that contains our global
459 // // node numbers and those of our neighbor.
460 // // If the set is the same size as the number
461 // // of nodes in both elements then they must
462 // // be connected to the same nodes.
463 // std::set<unsigned int> nodes_set;
464 
465 // for (unsigned int n=0; n<this->n_nodes(); n++)
466 // {
467 // nodes_set.insert(this->node(n));
468 // nodes_set.insert(rhs.node(n));
469 // }
470 
471 // // If this passes the elements are connected
472 // // to the same global nodes
473 // if (nodes_set.size() == this->n_nodes())
474 // return true;
475 // }
476 
477 // // If we get here it is because the elements either
478 // // do not have the same number of nodes or they are
479 // // connected to different nodes. Either way they
480 // // are not the same element
481 // return false;
482 
483  // Useful typedefs
484  typedef std::vector<dof_id_type>::iterator iterator;
485 
486 
487  // Elements can only be equal if they
488  // contain the same number of nodes.
489  // However, we will only test the vertices,
490  // which is sufficient & cheaper
491  if (this->n_nodes() == rhs_elem->n_nodes())
492  {
493  // The number of nodes in the element
494  const unsigned int nn = this->n_nodes();
495 
496  // Create a vector that contains our global
497  // node numbers and those of our neighbor.
498  // If the sorted, unique vector is the same size
499  // as the number of nodes in both elements then
500  // they must be connected to the same nodes.
501  //
502  // The vector will be no larger than 2*n_nodes(),
503  // so we might as well reserve the space.
504  std::vector<dof_id_type> common_nodes;
505  common_nodes.reserve (2*nn);
506 
507  // Add the global indices of the nodes
508  for (unsigned int n=0; n<nn; n++)
509  {
510  common_nodes.push_back (this->node(n));
511  common_nodes.push_back (rhs_elem->node(n));
512  }
513 
514  // Sort the vector and find out how long
515  // the sorted vector is.
516  std::sort (common_nodes.begin(), common_nodes.end());
517 
518  iterator new_end = std::unique (common_nodes.begin(),
519  common_nodes.end());
520 
521  const int new_size = libmesh_cast_int<int>
522  (std::distance (common_nodes.begin(), new_end));
523 
524  // If this passes the elements are connected
525  // to the same global vertex nodes
526  if (new_size == static_cast<int>(nn))
527  return true;
528  }
529 
530  // If we get here it is because the elements either
531  // do not have the same number of nodes or they are
532  // connected to different nodes. Either way they
533  // are not the same element
534  return false;
535 }
536 
537 
538 
539 bool Elem::is_semilocal(const processor_id_type my_pid) const
540 {
541  std::set<const Elem *> point_neighbors;
542 
543  this->find_point_neighbors(point_neighbors);
544 
545  std::set<const Elem*>::const_iterator it = point_neighbors.begin();
546  const std::set<const Elem*>::const_iterator end = point_neighbors.end();
547 
548  for (; it != end; ++it)
549  {
550  const Elem* elem = *it;
551  if (elem->processor_id() == my_pid)
552  return true;
553  }
554 
555  return false;
556 }
557 
558 
559 
560 bool Elem::contains_vertex_of(const Elem *e) const
561 {
562  // Our vertices are the first numbered nodes
563  for (unsigned int n = 0; n != e->n_vertices(); ++n)
564  if (this->contains_point(e->point(n)))
565  return true;
566  return false;
567 }
568 
569 
570 
571 bool Elem::contains_edge_of(const Elem *e) const
572 {
573  unsigned int num_contained_edges = 0;
574 
575  // Our vertices are the first numbered nodes
576  for (unsigned int n = 0; n != e->n_vertices(); ++n)
577  {
578  if (this->contains_point(e->point(n)))
579  {
580  num_contained_edges++;
581  if(num_contained_edges>=2)
582  {
583  return true;
584  }
585  }
586  }
587  return false;
588 }
589 
590 
591 
593  std::set<const Elem *> &neighbor_set) const
594 {
595  libmesh_assert(this->contains_point(p));
596 
597  neighbor_set.clear();
598  neighbor_set.insert(this);
599 
600  std::set<const Elem *> untested_set, next_untested_set;
601  untested_set.insert(this);
602 
603  while (!untested_set.empty())
604  {
605  // Loop over all the elements in the patch that haven't already
606  // been tested
607  std::set<const Elem*>::const_iterator it = untested_set.begin();
608  const std::set<const Elem*>::const_iterator end = untested_set.end();
609 
610  for (; it != end; ++it)
611  {
612  const Elem* elem = *it;
613 
614  for (unsigned int s=0; s<elem->n_sides(); s++)
615  {
616  const Elem* current_neighbor = elem->neighbor(s);
617  if (current_neighbor &&
618  current_neighbor != remote_elem) // we have a real neighbor on this side
619  {
620  if (current_neighbor->active()) // ... if it is active
621  {
622  if (current_neighbor->contains_point(p)) // ... and touches p
623  {
624  // Make sure we'll test it
625  if (!neighbor_set.count(current_neighbor))
626  next_untested_set.insert (current_neighbor);
627 
628  // And add it
629  neighbor_set.insert (current_neighbor);
630  }
631  }
632 #ifdef LIBMESH_ENABLE_AMR
633  else // ... the neighbor is *not* active,
634  { // ... so add *all* neighboring
635  // active children that touch p
636  std::vector<const Elem*> active_neighbor_children;
637 
638  current_neighbor->active_family_tree_by_neighbor
639  (active_neighbor_children, elem);
640 
641  std::vector<const Elem*>::const_iterator
642  child_it = active_neighbor_children.begin();
643  const std::vector<const Elem*>::const_iterator
644  child_end = active_neighbor_children.end();
645  for (; child_it != child_end; ++child_it)
646  {
647  const Elem *current_child = *child_it;
648  if (current_child->contains_point(p))
649  {
650  // Make sure we'll test it
651  if (!neighbor_set.count(current_child))
652  next_untested_set.insert (current_child);
653 
654  neighbor_set.insert (current_child);
655  }
656  }
657  }
658 #endif // #ifdef LIBMESH_ENABLE_AMR
659  }
660  }
661  }
662  untested_set.swap(next_untested_set);
663  next_untested_set.clear();
664  }
665 }
666 
667 
668 
669 void Elem::find_point_neighbors(std::set<const Elem *> &neighbor_set) const
670 {
671  neighbor_set.clear();
672  neighbor_set.insert(this);
673 
674  std::set<const Elem *> untested_set, next_untested_set;
675  untested_set.insert(this);
676 
677  while (!untested_set.empty())
678  {
679  // Loop over all the elements in the patch that haven't already
680  // been tested
681  std::set<const Elem*>::const_iterator it = untested_set.begin();
682  const std::set<const Elem*>::const_iterator end = untested_set.end();
683 
684  for (; it != end; ++it)
685  {
686  const Elem* elem = *it;
687 
688  for (unsigned int s=0; s<elem->n_sides(); s++)
689  {
690  const Elem* current_neighbor = elem->neighbor(s);
691  if (current_neighbor &&
692  current_neighbor != remote_elem) // we have a real neighbor on this side
693  {
694  if (current_neighbor->active()) // ... if it is active
695  {
696  if (this->contains_vertex_of(current_neighbor) // ... and touches us
697  || current_neighbor->contains_vertex_of(this))
698  {
699  // Make sure we'll test it
700  if (!neighbor_set.count(current_neighbor))
701  next_untested_set.insert (current_neighbor);
702 
703  // And add it
704  neighbor_set.insert (current_neighbor);
705  }
706  }
707 #ifdef LIBMESH_ENABLE_AMR
708  else // ... the neighbor is *not* active,
709  { // ... so add *all* neighboring
710  // active children
711  std::vector<const Elem*> active_neighbor_children;
712 
713  current_neighbor->active_family_tree_by_neighbor
714  (active_neighbor_children, elem);
715 
716  std::vector<const Elem*>::const_iterator
717  child_it = active_neighbor_children.begin();
718  const std::vector<const Elem*>::const_iterator
719  child_end = active_neighbor_children.end();
720  for (; child_it != child_end; ++child_it)
721  {
722  const Elem *current_child = *child_it;
723  if (this->contains_vertex_of(current_child) ||
724  (current_child)->contains_vertex_of(this))
725  {
726  // Make sure we'll test it
727  if (!neighbor_set.count(current_child))
728  next_untested_set.insert (current_child);
729 
730  neighbor_set.insert (current_child);
731  }
732  }
733  }
734 #endif // #ifdef LIBMESH_ENABLE_AMR
735  }
736  }
737  }
738  untested_set.swap(next_untested_set);
739  next_untested_set.clear();
740  }
741 }
742 
743 
744 
746  const Point& p2,
747  std::set<const Elem *> &neighbor_set) const
748 {
749  // Simple but perhaps suboptimal code: find elements containing the
750  // first point, then winnow this set down by removing elements which
751  // don't also contain the second point
752 
753  libmesh_assert(this->contains_point(p2));
754  this->find_point_neighbors(p1, neighbor_set);
755 
756  std::set<const Elem*>::iterator it = neighbor_set.begin();
757  const std::set<const Elem*>::iterator end = neighbor_set.end();
758 
759  while(it != end) {
760  std::set<const Elem*>::iterator current = it++;
761 
762  const Elem* elem = *current;
763  // This won't invalidate iterator it, because it is already
764  // pointing to the next element
765  if (!elem->contains_point(p2))
766  neighbor_set.erase(current);
767  }
768 }
769 
770 
771 
772 void Elem::find_edge_neighbors(std::set<const Elem *> &neighbor_set) const
773 {
774  neighbor_set.clear();
775  neighbor_set.insert(this);
776 
777  std::set<const Elem *> untested_set, next_untested_set;
778  untested_set.insert(this);
779 
780  while (!untested_set.empty())
781  {
782  // Loop over all the elements in the patch that haven't already
783  // been tested
784  std::set<const Elem*>::const_iterator it = untested_set.begin();
785  const std::set<const Elem*>::const_iterator end = untested_set.end();
786 
787  for (; it != end; ++it)
788  {
789  const Elem* elem = *it;
790 
791  for (unsigned int s=0; s<elem->n_sides(); s++)
792  {
793  const Elem* current_neighbor = elem->neighbor(s);
794  if (current_neighbor &&
795  current_neighbor != remote_elem) // we have a real neighbor on this side
796  {
797  if (current_neighbor->active()) // ... if it is active
798  {
799  if (this->contains_edge_of(current_neighbor) // ... and touches us
800  || current_neighbor->contains_edge_of(this))
801  {
802  // Make sure we'll test it
803  if (!neighbor_set.count(current_neighbor))
804  next_untested_set.insert (current_neighbor);
805 
806  // And add it
807  neighbor_set.insert (current_neighbor);
808  }
809  }
810 #ifdef LIBMESH_ENABLE_AMR
811  else // ... the neighbor is *not* active,
812  { // ... so add *all* neighboring
813  // active children
814  std::vector<const Elem*> active_neighbor_children;
815 
816  current_neighbor->active_family_tree_by_neighbor
817  (active_neighbor_children, elem);
818 
819  std::vector<const Elem*>::const_iterator
820  child_it = active_neighbor_children.begin();
821  const std::vector<const Elem*>::const_iterator
822  child_end = active_neighbor_children.end();
823  for (; child_it != child_end; ++child_it)
824  {
825  const Elem *current_child = *child_it;
826  if (this->contains_edge_of(*child_it) ||
827  (*child_it)->contains_edge_of(this))
828  {
829  // Make sure we'll test it
830  if (!neighbor_set.count(current_child))
831  next_untested_set.insert (current_child);
832 
833  neighbor_set.insert (current_child);
834  }
835  }
836  }
837 #endif // #ifdef LIBMESH_ENABLE_AMR
838  }
839  }
840  }
841  untested_set.swap(next_untested_set);
842  next_untested_set.clear();
843  }
844 }
845 
846 #ifdef LIBMESH_ENABLE_PERIODIC
847 
848 Elem* Elem::topological_neighbor (const unsigned int i,
849  MeshBase & mesh,
850  const PointLocatorBase& point_locator,
851  const PeriodicBoundaries * pb)
852 {
853  libmesh_assert_less (i, this->n_neighbors());
854 
855  Elem * neighbor_i = this->neighbor(i);
856  if (neighbor_i != NULL)
857  return neighbor_i;
858 
859  if (pb)
860  {
861  // Since the neighbor is NULL it must be on a boundary. We need
862  // see if this is a periodic boundary in which case it will have a
863  // topological neighbor
864 
865  std::vector<boundary_id_type> boundary_ids = mesh.boundary_info->boundary_ids(this, i);
866  for (std::vector<boundary_id_type>::iterator j = boundary_ids.begin(); j != boundary_ids.end(); ++j)
867  if (pb->boundary(*j))
868  {
869  // Since the point locator inside of periodic boundaries
870  // returns a const pointer we will retrieve the proper
871  // pointer directly from the mesh object. Also since coarse
872  // elements do not have more refined neighbors we need to make
873  // sure that we don't return one of these types of neighbors.
874  neighbor_i = mesh.elem(pb->neighbor(*j, point_locator, this, i)->id());
875  if (level() < neighbor_i->level())
876  neighbor_i = neighbor_i->parent();
877  return neighbor_i;
878  }
879  }
880 
881  return NULL;
882 }
883 
884 
885 
886 const Elem* Elem::topological_neighbor (const unsigned int i,
887  const MeshBase & mesh,
888  const PointLocatorBase& point_locator,
889  const PeriodicBoundaries * pb) const
890 {
891  libmesh_assert_less (i, this->n_neighbors());
892 
893  const Elem * neighbor_i = this->neighbor(i);
894  if (neighbor_i != NULL)
895  return neighbor_i;
896 
897  if (pb)
898  {
899  // Since the neighbor is NULL it must be on a boundary. We need
900  // see if this is a periodic boundary in which case it will have a
901  // topological neighbor
902 
903  std::vector<boundary_id_type> boundary_ids = mesh.boundary_info->boundary_ids(this, i);
904  for (std::vector<boundary_id_type>::iterator j = boundary_ids.begin(); j != boundary_ids.end(); ++j)
905  if (pb->boundary(*j))
906  {
907  // Since the point locator inside of periodic boundaries
908  // returns a const pointer we will retrieve the proper
909  // pointer directly from the mesh object. Also since coarse
910  // elements do not have more refined neighbors we need to make
911  // sure that we don't return one of these types of neighbors.
912  neighbor_i = mesh.elem(pb->neighbor(*j, point_locator, this, i)->id());
913  if (level() < neighbor_i->level())
914  neighbor_i = neighbor_i->parent();
915  return neighbor_i;
916  }
917  }
918 
919  return NULL;
920 }
921 
922 
924  const MeshBase & mesh,
925  const PointLocatorBase& point_locator,
926  PeriodicBoundaries * pb) const
927 {
928  // First see if this is a normal "interior" neighbor
929  if (has_neighbor(elem))
930  return true;
931 
932  for (unsigned int n=0; n<this->n_neighbors(); n++)
933  if (this->topological_neighbor(n, mesh, point_locator, pb))
934  return true;
935 
936  return false;
937 }
938 
939 
940 #endif
941 
942 #ifdef DEBUG
943 
945 {
946  libmesh_assert(this->valid_id());
947  for (unsigned int n=0; n != this->n_nodes(); ++n)
948  {
949  libmesh_assert(this->get_node(n));
950  libmesh_assert(this->get_node(n)->valid_id());
951  }
952 }
953 
954 
955 
957 {
958  for (unsigned int s=0; s<this->n_neighbors(); s++)
959  {
960  const Elem *neigh = this->neighbor(s);
961 
962  // Any element might have a remote neighbor; checking
963  // to make sure that's not inaccurate is tough.
964  if (neigh == remote_elem)
965  continue;
966 
967  if (neigh)
968  {
969  // Only subactive elements have subactive neighbors
970  libmesh_assert (this->subactive() || !neigh->subactive());
971 
972  const Elem *elem = this;
973 
974  // If we're subactive but our neighbor isn't, its
975  // return neighbor link will be to our first active
976  // ancestor OR to our inactive ancestor of the same
977  // level as neigh,
978  if (this->subactive() && !neigh->subactive())
979  {
980  for (elem = this; !elem->active();
981  elem = elem->parent())
982  libmesh_assert(elem);
983  }
984  else
985  {
986  unsigned int rev = neigh->which_neighbor_am_i(elem);
987  libmesh_assert_less (rev, neigh->n_neighbors());
988 
989  if (this->subactive() && !neigh->subactive())
990  {
991  while (neigh->neighbor(rev) != elem)
992  {
993  libmesh_assert(elem->parent());
994  elem = elem->parent();
995  }
996  }
997  else
998  {
999  Elem *nn = neigh->neighbor(rev);
1000  libmesh_assert(nn);
1001 
1002  for (; elem != nn; elem = elem->parent())
1003  libmesh_assert(elem);
1004  }
1005  }
1006  }
1007  // If we don't have a neighbor and we're not subactive, our
1008  // ancestors shouldn't have any neighbors in this same
1009  // direction.
1010  else if (!this->subactive())
1011  {
1012  const Elem *my_parent = this->parent();
1013  if (my_parent &&
1014  // A parent with a different dimension isn't really one of
1015  // our ancestors, it means we're on a boundary mesh and this
1016  // is an interior mesh element for which we're on a side.
1017  // Nothing to test for in that case.
1018  (my_parent->dim() == this->dim()))
1019  libmesh_assert (!my_parent->neighbor(s));
1020  }
1021  }
1022 }
1023 
1024 #endif // DEBUG
1025 
1026 
1027 
1028 void Elem::make_links_to_me_local(unsigned int n)
1029 {
1030  Elem *neigh = this->neighbor(n);
1031 
1032  // Don't bother calling this function unless it's necessary
1033  libmesh_assert(neigh);
1034  libmesh_assert(!neigh->is_remote());
1035 
1036  // We never have neighbors more refined than us
1037  libmesh_assert_less_equal (neigh->level(), this->level());
1038 
1039  // We never have subactive neighbors of non subactive elements
1040  libmesh_assert(!neigh->subactive() || this->subactive());
1041 
1042  // If we have a neighbor less refined than us then it must not
1043  // have any more refined active descendants we could have
1044  // pointed to instead.
1045  libmesh_assert(neigh->level() == this->level() ||
1046  neigh->active());
1047 
1048  // If neigh is at our level, then its family might have
1049  // remote_elem neighbor links which need to point to us
1050  // instead, but if not, then we're done.
1051  if (neigh->level() != this->level())
1052  return;
1053 
1054  // If neigh is subactive then we're not updating its neighbor links
1055  // FIXME - this needs to change when we start using subactive
1056  // elements for more than just the two-phase
1057  // restriction/prolongation projections.
1058  if (neigh->subactive())
1059  return;
1060 
1061  // What side of neigh are we on? We can't use the usual Elem
1062  // method because we're in the middle of restoring topology
1063  const AutoPtr<Elem> my_side = this->side(n);
1064  unsigned int nn = 0;
1065  for (; nn != neigh->n_sides(); ++nn)
1066  {
1067  const AutoPtr<Elem> neigh_side = neigh->side(nn);
1068  if (*my_side == *neigh_side)
1069  break;
1070  }
1071 
1072  // we had better be on *some* side of neigh
1073  libmesh_assert_less (nn, neigh->n_sides());
1074 
1075  // Find any elements that ought to point to elem
1076  std::vector<const Elem*> neigh_family;
1077 #ifdef LIBMESH_ENABLE_AMR
1078  if (this->active())
1079  neigh->family_tree_by_side(neigh_family, nn);
1080  else
1081 #endif
1082  neigh_family.push_back(neigh);
1083 
1084  // And point them to elem
1085  for (unsigned int i = 0; i != neigh_family.size(); ++i)
1086  {
1087  Elem* neigh_family_member = const_cast<Elem*>(neigh_family[i]);
1088 
1089  // Ideally, the neighbor link ought to either be correct
1090  // already or ought to be to remote_elem.
1091  //
1092  // However, if we're redistributing a newly created elem,
1093  // after an AMR step but before find_neighbors has fixed up
1094  // neighbor links, we might have an out of date neighbor
1095  // link to elem's parent instead.
1096 #ifdef LIBMESH_ENABLE_AMR
1097  libmesh_assert((neigh_family_member->neighbor(nn) == this) ||
1098  (neigh_family_member->neighbor(nn) == remote_elem)
1099  || ((this->refinement_flag() == JUST_REFINED) &&
1100  (this->parent() != NULL) &&
1101  (neigh_family_member->neighbor(nn) == this->parent())));
1102 #else
1103  libmesh_assert((neigh_family_member->neighbor(nn) == this) ||
1104  (neigh_family_member->neighbor(nn) == remote_elem));
1105 #endif
1106 
1107  neigh_family_member->set_neighbor(nn, this);
1108  }
1109 }
1110 
1111 
1113 {
1114  libmesh_assert_not_equal_to (this, remote_elem);
1115 
1116  // We need to have handled any children first
1117 #if defined(LIBMESH_ENABLE_AMR) && defined(DEBUG)
1118  if (this->has_children())
1119  for (unsigned int c = 0; c != this->n_children(); ++c)
1120  {
1121  Elem *current_child = this->child(c);
1122  libmesh_assert_equal_to (current_child, remote_elem);
1123  }
1124 #endif
1125 
1126  // Remotify any neighbor links to non-subactive elements
1127  if (!this->subactive())
1128  {
1129  for (unsigned int s = 0; s != this->n_sides(); ++s)
1130  {
1131  Elem *neigh = this->neighbor(s);
1132  if (neigh && neigh != remote_elem && !neigh->subactive())
1133  {
1134  // My neighbor should never be more refined than me; my real
1135  // neighbor would have been its parent in that case.
1136  libmesh_assert_greater_equal (this->level(), neigh->level());
1137 
1138  if (this->level() == neigh->level() &&
1139  neigh->has_neighbor(this))
1140  {
1141 #ifdef LIBMESH_ENABLE_AMR
1142  // My neighbor may have descendants which also consider me a
1143  // neighbor
1144  std::vector<const Elem*> family;
1145  neigh->family_tree_by_neighbor (family, this);
1146 
1147  // FIXME - There's a lot of ugly const_casts here; we
1148  // may want to make remote_elem non-const and create
1149  // non-const versions of the family_tree methods
1150  for (unsigned int i=0; i != family.size(); ++i)
1151  {
1152  Elem *n = const_cast<Elem*>(family[i]);
1153  libmesh_assert (n);
1154  if (n == remote_elem)
1155  continue;
1156  unsigned int my_s = n->which_neighbor_am_i(this);
1157  libmesh_assert_less (my_s, n->n_neighbors());
1158  libmesh_assert_equal_to (n->neighbor(my_s), this);
1159  n->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
1160  }
1161 #else
1162  unsigned int my_s = neigh->which_neighbor_am_i(this);
1163  libmesh_assert_less (my_s, neigh->n_neighbors());
1164  libmesh_assert_equal_to (neigh->neighbor(my_s), this);
1165  neigh->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
1166 #endif
1167  }
1168 #ifdef LIBMESH_ENABLE_AMR
1169  // Even if my neighbor doesn't link back to me, it might
1170  // have subactive descendants which do
1171  else if (neigh->has_children())
1172  {
1173  // If my neighbor at the same level doesn't have me as a
1174  // neighbor, I must be subactive
1175  libmesh_assert(this->level() > neigh->level() ||
1176  this->subactive());
1177 
1178  // My neighbor must have some ancestor of mine as a
1179  // neighbor
1180  Elem *my_ancestor = this->parent();
1181  libmesh_assert(my_ancestor);
1182  while (!neigh->has_neighbor(my_ancestor))
1183  {
1184  my_ancestor = my_ancestor->parent();
1185  libmesh_assert(my_ancestor);
1186  }
1187 
1188  // My neighbor may have descendants which consider me a
1189  // neighbor
1190  std::vector<const Elem*> family;
1191  neigh->family_tree_by_subneighbor (family, my_ancestor, this);
1192 
1193  // FIXME - There's a lot of ugly const_casts here; we
1194  // may want to make remote_elem non-const and create
1195  // non-const versions of the family_tree methods
1196  for (unsigned int i=0; i != family.size(); ++i)
1197  {
1198  Elem *n = const_cast<Elem*>(family[i]);
1199  libmesh_assert (n);
1200  if (n == remote_elem)
1201  continue;
1202  unsigned int my_s = n->which_neighbor_am_i(this);
1203  libmesh_assert_less (my_s, n->n_neighbors());
1204  libmesh_assert_equal_to (n->neighbor(my_s), this);
1205  n->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
1206  }
1207  }
1208 #endif
1209  }
1210  }
1211  }
1212 
1213 #ifdef LIBMESH_ENABLE_AMR
1214  // Remotify parent's child link
1215  Elem *my_parent = this->parent();
1216  if (my_parent &&
1217  // As long as it's not already remote
1218  my_parent != remote_elem &&
1219  // And it's a real parent, not an interior parent
1220  this->dim() == my_parent->dim())
1221  {
1222  unsigned int me = my_parent->which_child_am_i(this);
1223  libmesh_assert_equal_to (my_parent->child(me), this);
1224  my_parent->set_child(me, const_cast<RemoteElem*>(remote_elem));
1225  }
1226 #endif
1227 }
1228 
1229 
1230 
1231 void Elem::write_connectivity (std::ostream& out_stream,
1232  const IOPackage iop) const
1233 {
1234  libmesh_assert (out_stream.good());
1236  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
1237 
1238  switch (iop)
1239  {
1240  case TECPLOT:
1241  {
1242  // This connectivity vector will be used repeatedly instead
1243  // of being reconstructed inside the loop.
1244  std::vector<dof_id_type> conn;
1245  for (unsigned int sc=0; sc <this->n_sub_elem(); sc++)
1246  {
1247  this->connectivity(sc, TECPLOT, conn);
1248 
1249  std::copy(conn.begin(),
1250  conn.end(),
1251  std::ostream_iterator<dof_id_type>(out_stream, " "));
1252 
1253  out_stream << '\n';
1254  }
1255  return;
1256  }
1257 
1258  case UCD:
1259  {
1260  for (unsigned int i=0; i<this->n_nodes(); i++)
1261  out_stream << this->node(i)+1 << "\t";
1262 
1263  out_stream << '\n';
1264  return;
1265  }
1266 
1267  default:
1268  libmesh_error();
1269  }
1270 
1271  libmesh_error();
1272 }
1273 
1274 
1275 // void Elem::write_tecplot_connectivity(std::ostream& out) const
1276 // {
1277 // libmesh_assert (!out.bad());
1278 // libmesh_assert(_nodes);
1279 
1280 // // This connectivity vector will be used repeatedly instead
1281 // // of being reconstructed inside the loop.
1282 // std::vector<dof_id_type> conn;
1283 // for (unsigned int sc=0; sc <this->n_sub_elem(); sc++)
1284 // {
1285 // this->connectivity(sc, TECPLOT, conn);
1286 
1287 // std::copy(conn.begin(),
1288 // conn.end(),
1289 // std::ostream_iterator<dof_id_type>(out, " "));
1290 
1291 // out << std::endl;
1292 // }
1293 // }
1294 
1295 
1296 
1297 // void Elem::write_ucd_connectivity(std::ostream &out) const
1298 // {
1299 // libmesh_assert (out);
1300 // libmesh_assert(_nodes);
1301 
1302 // for (unsigned int i=0; i<this->n_nodes(); i++)
1303 // out << this->node(i)+1 << "\t";
1304 
1305 // out << std::endl;
1306 // }
1307 
1308 
1309 
1311 {
1312  switch (q)
1313  {
1317  default:
1318  {
1319  libmesh_do_once( libmesh_here();
1320 
1321  libMesh::err << "ERROR: unknown quality metric: "
1323  << std::endl
1324  << "Cowardly returning 1."
1325  << std::endl; );
1326 
1327  return 1.;
1328  }
1329  }
1330 
1331 
1332  // Will never get here...
1333  libmesh_error();
1334  return 0.;
1335 }
1336 
1337 
1338 
1339 bool Elem::ancestor() const
1340 {
1341 #ifdef LIBMESH_ENABLE_AMR
1342 
1343  if (this->active())
1344  return false;
1345 
1346 if (!this->has_children())
1347  return false;
1348  if (this->child(0)->active())
1349  return true;
1350 
1351  return this->child(0)->ancestor();
1352 #else
1353  return false;
1354 #endif
1355 }
1356 
1357 
1358 
1359 #ifdef LIBMESH_ENABLE_AMR
1360 
1361 void Elem::add_child (Elem* elem)
1362 {
1363  if(_children == NULL)
1364  {
1365  _children = new Elem*[this->n_children()];
1366 
1367  for (unsigned int c=0; c<this->n_children(); c++)
1368  this->set_child(c, NULL);
1369  }
1370 
1371  for (unsigned int c=0; c<this->n_children(); c++)
1372  {
1373  if(this->_children[c] == NULL || this->_children[c] == remote_elem)
1374  {
1375  libmesh_assert_equal_to (this, elem->parent());
1376  this->set_child(c, elem);
1377  return;
1378  }
1379  }
1380 
1381  libMesh::err << "Error: Tried to add a child to an element with full children array"
1382  << std::endl;
1383  libmesh_error();
1384 }
1385 
1386 
1387 
1388 void Elem::add_child (Elem* elem, unsigned int c)
1389 {
1390  if(!this->has_children())
1391  {
1392  _children = new Elem*[this->n_children()];
1393 
1394  for (unsigned int i=0; i<this->n_children(); i++)
1395  this->set_child(i, NULL);
1396  }
1397 
1398  libmesh_assert (this->_children[c] == NULL || this->child(c) == remote_elem);
1399  libmesh_assert (elem == remote_elem || this == elem->parent());
1400 
1401  this->set_child(c, elem);
1402 }
1403 
1404 
1405 
1406 void Elem::replace_child (Elem* elem, unsigned int c)
1407 {
1408  libmesh_assert(this->has_children());
1409 
1410  libmesh_assert(this->child(c));
1411 
1412  this->set_child(c, elem);
1413 }
1414 
1415 
1416 
1417 bool Elem::is_child_on_edge(const unsigned int libmesh_dbg_var(c),
1418  const unsigned int e) const
1419 {
1420  libmesh_assert_less (c, this->n_children());
1421  libmesh_assert_less (e, this->n_edges());
1422 
1423  AutoPtr<Elem> my_edge = this->build_edge(e);
1424  AutoPtr<Elem> child_edge = this->build_edge(e);
1425 
1426  // We're assuming that an overlapping child edge has the same
1427  // number and orientation as its parent
1428  return (child_edge->node(0) == my_edge->node(0) ||
1429  child_edge->node(1) == my_edge->node(1));
1430 }
1431 
1432 
1433 void Elem::family_tree (std::vector<const Elem*>& family,
1434  const bool reset) const
1435 {
1436  // The "family tree" doesn't include subactive elements
1437  libmesh_assert(!this->subactive());
1438 
1439  // Clear the vector if the flag reset tells us to.
1440  if (reset)
1441  family.clear();
1442 
1443  // Add this element to the family tree.
1444  family.push_back(this);
1445 
1446  // Recurse into the elements children, if it has them.
1447  // Do not clear the vector any more.
1448  if (!this->active())
1449  for (unsigned int c=0; c<this->n_children(); c++)
1450  if (!this->child(c)->is_remote())
1451  this->child(c)->family_tree (family, false);
1452 }
1453 
1454 
1455 
1456 void Elem::total_family_tree (std::vector<const Elem*>& family,
1457  const bool reset) const
1458 {
1459  // Clear the vector if the flag reset tells us to.
1460  if (reset)
1461  family.clear();
1462 
1463  // Add this element to the family tree.
1464  family.push_back(this);
1465 
1466  // Recurse into the elements children, if it has them.
1467  // Do not clear the vector any more.
1468  if (this->has_children())
1469  for (unsigned int c=0; c<this->n_children(); c++)
1470  if (!this->child(c)->is_remote())
1471  this->child(c)->total_family_tree (family, false);
1472 }
1473 
1474 
1475 
1476 void Elem::active_family_tree (std::vector<const Elem*>& active_family,
1477  const bool reset) const
1478 {
1479  // The "family tree" doesn't include subactive elements
1480  libmesh_assert(!this->subactive());
1481 
1482  // Clear the vector if the flag reset tells us to.
1483  if (reset)
1484  active_family.clear();
1485 
1486  // Add this element to the family tree if it is active
1487  if (this->active())
1488  active_family.push_back(this);
1489 
1490  // Otherwise recurse into the element's children.
1491  // Do not clear the vector any more.
1492  else
1493  for (unsigned int c=0; c<this->n_children(); c++)
1494  if (!this->child(c)->is_remote())
1495  this->child(c)->active_family_tree (active_family, false);
1496 }
1497 
1498 
1499 
1500 void Elem::family_tree_by_side (std::vector<const Elem*>& family,
1501  const unsigned int s,
1502  const bool reset) const
1503 {
1504  // The "family tree" doesn't include subactive elements
1505  libmesh_assert(!this->subactive());
1506 
1507  // Clear the vector if the flag reset tells us to.
1508  if (reset)
1509  family.clear();
1510 
1511  libmesh_assert_less (s, this->n_sides());
1512 
1513  // Add this element to the family tree.
1514  family.push_back(this);
1515 
1516  // Recurse into the elements children, if it has them.
1517  // Do not clear the vector any more.
1518  if (!this->active())
1519  for (unsigned int c=0; c<this->n_children(); c++)
1520  if (!this->child(c)->is_remote() && this->is_child_on_side(c, s))
1521  this->child(c)->family_tree_by_side (family, s, false);
1522 }
1523 
1524 
1525 
1526 void Elem::active_family_tree_by_side (std::vector<const Elem*>& family,
1527  const unsigned int s,
1528  const bool reset) const
1529 {
1530  // The "family tree" doesn't include subactive elements
1531  libmesh_assert(!this->subactive());
1532 
1533  // Clear the vector if the flag reset tells us to.
1534  if (reset)
1535  family.clear();
1536 
1537  libmesh_assert_less (s, this->n_sides());
1538 
1539  // Add an active element to the family tree.
1540  if (this->active())
1541  family.push_back(this);
1542 
1543  // Or recurse into an ancestor element's children.
1544  // Do not clear the vector any more.
1545  else
1546  for (unsigned int c=0; c<this->n_children(); c++)
1547  if (!this->child(c)->is_remote() && this->is_child_on_side(c, s))
1548  this->child(c)->active_family_tree_by_side (family, s, false);
1549 }
1550 
1551 
1552 
1553 void Elem::family_tree_by_neighbor (std::vector<const Elem*>& family,
1554  const Elem* neighbor_in,
1555  const bool reset) const
1556 {
1557  // The "family tree" doesn't include subactive elements
1558  libmesh_assert(!this->subactive());
1559 
1560  // Clear the vector if the flag reset tells us to.
1561  if (reset)
1562  family.clear();
1563 
1564  // This only makes sense if we're already a neighbor
1565  libmesh_assert (this->has_neighbor(neighbor_in));
1566 
1567  // Add this element to the family tree.
1568  family.push_back(this);
1569 
1570  // Recurse into the elements children, if it's not active.
1571  // Do not clear the vector any more.
1572  if (!this->active())
1573  for (unsigned int c=0; c<this->n_children(); c++)
1574  {
1575  Elem *current_child = this->child(c);
1576  if (current_child != remote_elem && current_child->has_neighbor(neighbor_in))
1577  current_child->family_tree_by_neighbor (family, neighbor_in, false);
1578  }
1579 }
1580 
1581 
1582 
1583 void Elem::family_tree_by_subneighbor (std::vector<const Elem*>& family,
1584  const Elem* neighbor_in,
1585  const Elem* subneighbor,
1586  const bool reset) const
1587 {
1588  // The "family tree" doesn't include subactive elements
1589  libmesh_assert(!this->subactive());
1590 
1591  // Clear the vector if the flag reset tells us to.
1592  if (reset)
1593  family.clear();
1594 
1595  // To simplifly this function we need an existing neighbor
1596  libmesh_assert (neighbor_in);
1597  libmesh_assert_not_equal_to (neighbor_in, remote_elem);
1598  libmesh_assert (this->has_neighbor(neighbor_in));
1599 
1600  // This only makes sense if subneighbor descends from neighbor
1601  libmesh_assert (subneighbor);
1602  libmesh_assert_not_equal_to (subneighbor, remote_elem);
1603  libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
1604 
1605  // Add this element to the family tree if applicable.
1606  if (neighbor_in == subneighbor)
1607  family.push_back(this);
1608 
1609  // Recurse into the elements children, if it's not active.
1610  // Do not clear the vector any more.
1611  if (!this->active())
1612  for (unsigned int c=0; c != this->n_children(); ++c)
1613  {
1614  Elem *current_child = this->child(c);
1615  if (current_child != remote_elem)
1616  for (unsigned int s=0; s != current_child->n_sides(); ++s)
1617  {
1618  Elem *child_neigh = current_child->neighbor(s);
1619  if (child_neigh &&
1620  (child_neigh == neighbor_in ||
1621  (child_neigh->parent() == neighbor_in &&
1622  child_neigh->is_ancestor_of(subneighbor))))
1623  current_child->family_tree_by_subneighbor (family, child_neigh,
1624  subneighbor, false);
1625  }
1626  }
1627 }
1628 
1629 
1630 
1631 void Elem::active_family_tree_by_neighbor (std::vector<const Elem*>& family,
1632  const Elem* neighbor_in,
1633  const bool reset) const
1634 {
1635  // The "family tree" doesn't include subactive elements
1636  libmesh_assert(!this->subactive());
1637 
1638  // Clear the vector if the flag reset tells us to.
1639  if (reset)
1640  family.clear();
1641 
1642  // This only makes sense if we're already a neighbor
1643  if (this->level() >= neighbor_in->level())
1644  libmesh_assert (this->has_neighbor(neighbor_in));
1645 
1646  // Add an active element to the family tree.
1647  if (this->active())
1648  family.push_back(this);
1649 
1650  // Or recurse into an ancestor element's children.
1651  // Do not clear the vector any more.
1652  else if (!this->active())
1653  for (unsigned int c=0; c<this->n_children(); c++)
1654  {
1655  Elem *current_child = this->child(c);
1656  if (current_child != remote_elem && current_child->has_neighbor(neighbor_in))
1657  current_child->active_family_tree_by_neighbor (family, neighbor_in, false);
1658  }
1659 }
1660 
1661 
1662 
1663 unsigned int Elem::min_p_level_by_neighbor(const Elem* neighbor_in,
1664  unsigned int current_min) const
1665 {
1666  libmesh_assert(!this->subactive());
1667  libmesh_assert(neighbor_in->active());
1668 
1669  // If we're an active element this is simple
1670  if (this->active())
1671  return std::min(current_min, this->p_level());
1672 
1673  libmesh_assert(has_neighbor(neighbor_in));
1674 
1675  // The p_level() of an ancestor element is already the minimum
1676  // p_level() of its children - so if that's high enough, we don't
1677  // need to examine any children.
1678  if (current_min <= this->p_level())
1679  return current_min;
1680 
1681  unsigned int min_p_level = current_min;
1682 
1683  for (unsigned int c=0; c<this->n_children(); c++)
1684  {
1685  const Elem* const current_child = this->child(c);
1686  if (current_child != remote_elem && current_child->has_neighbor(neighbor_in))
1687  min_p_level =
1688  current_child->min_p_level_by_neighbor(neighbor_in,
1689  min_p_level);
1690  }
1691 
1692  return min_p_level;
1693 }
1694 
1695 
1696 unsigned int Elem::min_new_p_level_by_neighbor(const Elem* neighbor_in,
1697  unsigned int current_min) const
1698 {
1699  libmesh_assert(!this->subactive());
1700  libmesh_assert(neighbor_in->active());
1701 
1702  // If we're an active element this is simple
1703  if (this->active())
1704  {
1705  unsigned int new_p_level = this->p_level();
1706  if (this->p_refinement_flag() == Elem::REFINE)
1707  new_p_level += 1;
1708  if (this->p_refinement_flag() == Elem::COARSEN)
1709  {
1710  libmesh_assert_greater (new_p_level, 0);
1711  new_p_level -= 1;
1712  }
1713  return std::min(current_min, new_p_level);
1714  }
1715 
1716  libmesh_assert(has_neighbor(neighbor_in));
1717 
1718  unsigned int min_p_level = current_min;
1719 
1720  for (unsigned int c=0; c<this->n_children(); c++)
1721  {
1722  const Elem* const current_child = this->child(c);
1723  if (current_child && current_child != remote_elem)
1724  if (current_child->has_neighbor(neighbor_in))
1725  min_p_level =
1726  current_child->min_new_p_level_by_neighbor(neighbor_in,
1727  min_p_level);
1728  }
1729 
1730  return min_p_level;
1731 }
1732 
1733 #endif // #ifdef LIBMESH_ENABLE_AMR
1734 
1735 
1736 
1737 
1738 bool Elem::contains_point (const Point& p, Real tol) const
1739 {
1740  // We currently allow the user to enlarge the bounding box by
1741  // providing a tol > TOLERANCE (so this routine is identical to
1742  // Elem::close_to_point()), but print a warning so that the
1743  // user can eventually switch his code over to calling close_to_point()
1744  // instead, which is intended to be used for this purpose.
1745  if ( tol > TOLERANCE )
1746  {
1747  libmesh_do_once(libMesh::err
1748  << "WARNING: Resizing bounding box to match user-specified tolerance!\n"
1749  << "In the future, calls to Elem::contains_point() with tol > TOLERANCE\n"
1750  << "will be more optimized, but should not be used\n"
1751  << "to search for points 'close to' elements!\n"
1752  << "Instead, use Elem::close_to_point() for this purpose.\n"
1753  << std::endl;);
1754  return this->point_test(p, tol, tol);
1755  }
1756  else
1757  return this->point_test(p, TOLERANCE, tol);
1758 }
1759 
1760 
1761 
1762 
1763 bool Elem::close_to_point (const Point& p, Real tol) const
1764 {
1765  // This test uses the user's passed-in tolerance for the
1766  // bounding box test as well, thereby allowing the routine to
1767  // find points which are not only "in" the element, but also
1768  // "nearby" to within some tolerance.
1769  return this->point_test(p, tol, tol);
1770 }
1771 
1772 
1773 
1774 
1775 bool Elem::point_test(const Point& p, Real box_tol, Real map_tol) const
1776 {
1777  libmesh_assert_greater (box_tol, 0.);
1778  libmesh_assert_greater (map_tol, 0.);
1779 
1780  // This is a great optimization on first order elements, but it
1781  // could return false negatives on higher orders
1782  if (this->default_order() == FIRST)
1783  {
1784  // Check to make sure the element *could* contain this point, so we
1785  // can avoid an expensive inverse_map call if it doesn't.
1786  bool
1787 #if LIBMESH_DIM > 2
1788  point_above_min_z = false,
1789  point_below_max_z = false,
1790 #endif
1791 #if LIBMESH_DIM > 1
1792  point_above_min_y = false,
1793  point_below_max_y = false,
1794 #endif
1795  point_above_min_x = false,
1796  point_below_max_x = false;
1797 
1798  // For relative bounding box checks in physical space
1799  const Real my_hmax = this->hmax();
1800 
1801  for (unsigned int n=0; n != this->n_nodes(); ++n)
1802  {
1803  Point pe = this->point(n);
1804  point_above_min_x = point_above_min_x || (pe(0) - my_hmax*box_tol <= p(0));
1805  point_below_max_x = point_below_max_x || (pe(0) + my_hmax*box_tol >= p(0));
1806 #if LIBMESH_DIM > 1
1807  point_above_min_y = point_above_min_y || (pe(1) - my_hmax*box_tol <= p(1));
1808  point_below_max_y = point_below_max_y || (pe(1) + my_hmax*box_tol >= p(1));
1809 #endif
1810 #if LIBMESH_DIM > 2
1811  point_above_min_z = point_above_min_z || (pe(2) - my_hmax*box_tol <= p(2));
1812  point_below_max_z = point_below_max_z || (pe(2) + my_hmax*box_tol >= p(2));
1813 #endif
1814  }
1815 
1816  if (
1817 #if LIBMESH_DIM > 2
1818  !point_above_min_z ||
1819  !point_below_max_z ||
1820 #endif
1821 #if LIBMESH_DIM > 1
1822  !point_above_min_y ||
1823  !point_below_max_y ||
1824 #endif
1825  !point_above_min_x ||
1826  !point_below_max_x)
1827  return false;
1828  }
1829 
1830  // Declare a basic FEType. Will be a Lagrange
1831  // element by default.
1832  FEType fe_type(this->default_order());
1833 
1834  // To be on the safe side, we converge the inverse_map() iteration
1835  // to a slightly tighter tolerance than that requested by the
1836  // user...
1837  const Point mapped_point = FEInterface::inverse_map(this->dim(),
1838  fe_type,
1839  this,
1840  p,
1841  0.1*map_tol, // <- this is |dx| tolerance, the Newton residual should be ~ |dx|^2
1842  /*secure=*/ false);
1843 
1844  // Check that the refspace point maps back to p! This is only necessary
1845  // for 1D and 2D elements, 3D elements always live in 3D.
1846  //
1847  // TODO: The contains_point() function could most likely be implemented
1848  // more efficiently in the element sub-classes themselves, at least for
1849  // the linear element types.
1850  if (this->dim() < 3)
1851  {
1852  Point xyz = FEInterface::map(this->dim(),
1853  fe_type,
1854  this,
1855  mapped_point);
1856 
1857  // Compute the distance between the original point and the re-mapped point.
1858  // They should be in the same place.
1859  Real dist = (xyz - p).size();
1860 
1861 
1862  // If dist is larger than some fraction of the tolerance, then return false.
1863  // This can happen when e.g. a 2D element is living in 3D, and
1864  // FEInterface::inverse_map() maps p onto the projection of the element,
1865  // effectively "tricking" FEInterface::on_reference_element().
1866  if (dist > this->hmax() * map_tol)
1867  return false;
1868  }
1869 
1870 
1871 
1872  return FEInterface::on_reference_element(mapped_point, this->type(), map_tol);
1873 }
1874 
1875 
1876 
1877 
1878 void Elem::print_info (std::ostream& os) const
1879 {
1880  os << this->get_info()
1881  << std::endl;
1882 }
1883 
1884 
1885 
1886 std::string Elem::get_info () const
1887 {
1888  std::ostringstream oss;
1889 
1890  oss << " Elem Information" << '\n'
1891  << " id()=";
1892 
1893  if (this->valid_id())
1894  oss << this->id();
1895  else
1896  oss << "invalid";
1897 
1898  oss << ", processor_id()=" << this->processor_id() << '\n';
1899 
1900  oss << " type()=" << Utility::enum_to_string(this->type()) << '\n'
1901  << " dim()=" << this->dim() << '\n'
1902  << " n_nodes()=" << this->n_nodes() << '\n';
1903 
1904  for (unsigned int n=0; n != this->n_nodes(); ++n)
1905  oss << " " << n << *this->get_node(n);
1906 
1907  oss << " n_sides()=" << this->n_sides() << '\n';
1908 
1909  for (unsigned int s=0; s != this->n_sides(); ++s)
1910  {
1911  oss << " neighbor(" << s << ")=";
1912  if (this->neighbor(s))
1913  oss << this->neighbor(s)->id() << '\n';
1914  else
1915  oss << "NULL\n";
1916  }
1917 
1918  oss << " hmin()=" << this->hmin()
1919  << ", hmax()=" << this->hmax() << '\n'
1920  << " volume()=" << this->volume() << '\n'
1921  << " active()=" << this->active()
1922  << ", ancestor()=" << this->ancestor()
1923  << ", subactive()=" << this->subactive()
1924  << ", has_children()=" << this->has_children() << '\n'
1925  << " parent()=";
1926  if (this->parent())
1927  oss << this->parent()->id() << '\n';
1928  else
1929  oss << "NULL\n";
1930  oss << " level()=" << this->level()
1931  << ", p_level()=" << this->p_level() << '\n'
1932 #ifdef LIBMESH_ENABLE_AMR
1933  << " refinement_flag()=" << Utility::enum_to_string(this->refinement_flag()) << '\n'
1934  << " p_refinement_flag()=" << Utility::enum_to_string(this->p_refinement_flag()) << '\n'
1935 #endif
1936 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1937  << " infinite()=" << this->infinite() << '\n';
1938  if (this->infinite())
1939  oss << " origin()=" << this->origin() << '\n'
1940 #endif
1941  ;
1942 
1943  oss << " DoFs=";
1944  for (unsigned int s=0; s != this->n_systems(); ++s)
1945  for (unsigned int v=0; v != this->n_vars(s); ++v)
1946  for (unsigned int c=0; c != this->n_comp(s,v); ++c)
1947  oss << '(' << s << '/' << v << '/' << this->dof_number(s,v,c) << ") ";
1948 
1949 
1950  return oss.str();
1951 }
1952 
1953 
1954 
1956 {
1957  // Tell any of my neighbors about my death...
1958  // Looks strange, huh?
1959  for (unsigned int n=0; n<this->n_neighbors(); n++)
1960  {
1961  Elem* current_neighbor = this->neighbor(n);
1962  if (current_neighbor && current_neighbor != remote_elem)
1963  {
1964  // Note: it is possible that I see the neighbor
1965  // (which is coarser than me)
1966  // but they don't see me, so avoid that case.
1967  if (current_neighbor->level() == this->level())
1968  {
1969  const unsigned int w_n_a_i = current_neighbor->which_neighbor_am_i(this);
1970  libmesh_assert_less (w_n_a_i, current_neighbor->n_neighbors());
1971  current_neighbor->set_neighbor(w_n_a_i, NULL);
1972  this->set_neighbor(n, NULL);
1973  }
1974  }
1975  }
1976 }
1977 
1978 
1979 
1980 
1981 unsigned int Elem::n_second_order_adjacent_vertices (const unsigned int) const
1982 {
1983  // for linear elements, always return 0
1984  return 0;
1985 }
1986 
1987 
1988 
1989 unsigned short int Elem::second_order_adjacent_vertex (const unsigned int,
1990  const unsigned int) const
1991 {
1992  // for linear elements, always return 0
1993  return 0;
1994 }
1995 
1996 
1997 
1998 std::pair<unsigned short int, unsigned short int>
1999 Elem::second_order_child_vertex (const unsigned int) const
2000 {
2001  // for linear elements, always return 0
2002  return std::pair<unsigned short int, unsigned short int>(0,0);
2003 }
2004 
2005 
2006 
2008 {
2009  switch (et)
2010  {
2011  case EDGE2:
2012  case EDGE3:
2013  case EDGE4:
2014  return EDGE2;
2015  case TRI3:
2016  case TRI6:
2017  return TRI3;
2018  case QUAD4:
2019  case QUAD8:
2020  case QUAD9:
2021  return QUAD4;
2022  case TET4:
2023  case TET10:
2024  return TET4;
2025  case HEX8:
2026  case HEX27:
2027  case HEX20:
2028  return HEX8;
2029  case PRISM6:
2030  case PRISM15:
2031  case PRISM18:
2032  return PRISM6;
2033  case PYRAMID5:
2034  case PYRAMID14:
2035  return PYRAMID5;
2036 
2037 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2038 
2039  case INFQUAD4:
2040  case INFQUAD6:
2041  return INFQUAD4;
2042  case INFHEX8:
2043  case INFHEX16:
2044  case INFHEX18:
2045  return INFHEX8;
2046  case INFPRISM6:
2047  case INFPRISM12:
2048  return INFPRISM6;
2049 
2050 #endif
2051 
2052  default:
2053  // unknown element
2054  return INVALID_ELEM;
2055  }
2056 }
2057 
2058 
2059 
2061  const bool full_ordered)
2062 {
2063  /* for second-order elements, always return \p INVALID_ELEM
2064  * since second-order elements should not be converted
2065  * into something else. Only linear elements should
2066  * return something sensible here
2067  */
2068  switch (et)
2069  {
2070  case EDGE2:
2071  {
2072  // full_ordered not relevant
2073  return EDGE3;
2074  }
2075 
2076  case TRI3:
2077  {
2078  // full_ordered not relevant
2079  return TRI6;
2080  }
2081 
2082  case QUAD4:
2083  {
2084  if (full_ordered)
2085  return QUAD9;
2086  else
2087  return QUAD8;
2088  }
2089 
2090  case TET4:
2091  {
2092  // full_ordered not relevant
2093  return TET10;
2094  }
2095 
2096  case HEX8:
2097  {
2098  // see below how this correlates with INFHEX8
2099  if (full_ordered)
2100  return HEX27;
2101  else
2102  return HEX20;
2103  }
2104 
2105  case PRISM6:
2106  {
2107  if (full_ordered)
2108  return PRISM18;
2109  else
2110  return PRISM15;
2111  }
2112 
2113  case PYRAMID5:
2114  {
2115  if (full_ordered)
2116  return PYRAMID14;
2117  else // PYRAMID13 to be added...
2118  libmesh_error();
2119 
2120  return INVALID_ELEM;
2121  }
2122 
2123 
2124 
2125 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2126 
2127  // infinite elements
2128  case INFEDGE2:
2129  {
2130  // libmesh_error();
2131  return INVALID_ELEM;
2132  }
2133 
2134  case INFQUAD4:
2135  {
2136  // full_ordered not relevant
2137  return INFQUAD6;
2138  }
2139 
2140  case INFHEX8:
2141  {
2142  /*
2143  * Note that this matches with \p Hex8:
2144  * For full-ordered, \p InfHex18 and \p Hex27
2145  * belong together, and for not full-ordered,
2146  * \p InfHex16 and \p Hex20 belong together.
2147  */
2148  if (full_ordered)
2149  return INFHEX18;
2150  else
2151  return INFHEX16;
2152  }
2153 
2154  case INFPRISM6:
2155  {
2156  // full_ordered not relevant
2157  return INFPRISM12;
2158  }
2159 
2160 #endif
2161 
2162 
2163  default:
2164  {
2165  // second-order element
2166  return INVALID_ELEM;
2167  }
2168  }
2169 }
2170 
2171 
2172 
2174 {
2176  return side_iterator(this->_first_side(), this->_last_side(), bsp);
2177 }
2178 
2179 
2180 
2181 
2183 {
2185  return side_iterator(this->_last_side(), this->_last_side(), bsp);
2186 }
2187 
2188 
2189 
2190 
2192 {
2193  // The default implementation builds a finite element of the correct
2194  // order and sums up the JxW contributions. This can be expensive,
2195  // so the various element types can overload this method and compute
2196  // the volume more efficiently.
2197  FEType fe_type (this->default_order() , LAGRANGE);
2198 
2199  AutoPtr<FEBase> fe (FEBase::build(this->dim(),
2200  fe_type));
2201 
2202  const std::vector<Real>& JxW = fe->get_JxW();
2203 
2204  // The default quadrature rule should integrate the mass matrix,
2205  // thus it should be plenty to compute the area
2206  QGauss qrule (this->dim(), fe_type.default_quadrature_order());
2207 
2208  fe->attach_quadrature_rule(&qrule);
2209 
2210  fe->reinit(this);
2211 
2212  Real vol=0.;
2213  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
2214  vol += JxW[qp];
2215 
2216  return vol;
2217 
2218 }
2219 
2220 
2221 // ------------------------------------------------------------
2222 // Elem::PackedElem static data
2223 const unsigned int Elem::PackedElem::header_size = 10;
2224 
2225 
2226 // Elem::PackedElem member funcions
2227 void Elem::PackedElem::pack (std::vector<largest_id_type> &conn, const Elem* elem)
2228 {
2229  libmesh_assert(elem);
2230 
2231  // we can do at least this good. note that hopefully in general
2232  // the user will already have reserved the full space, which will render
2233  // this redundant
2234  conn.reserve (conn.size() + elem->packed_size());
2235 
2236 #ifdef LIBMESH_ENABLE_AMR
2237  conn.push_back (static_cast<largest_id_type>(elem->level()));
2238  conn.push_back (static_cast<largest_id_type>(elem->p_level()));
2239  conn.push_back (static_cast<largest_id_type>(elem->refinement_flag()));
2240  conn.push_back (static_cast<largest_id_type>(elem->p_refinement_flag()));
2241 #else
2242  conn.push_back (0);
2243  conn.push_back (0);
2244  conn.push_back (0);
2245  conn.push_back (0);
2246 #endif
2247  conn.push_back (static_cast<largest_id_type>(elem->type()));
2248  conn.push_back (elem->processor_id());
2249  conn.push_back (elem->subdomain_id());
2250  conn.push_back (elem->id());
2251 
2252 #ifdef LIBMESH_ENABLE_AMR
2253  // use parent_ID of -1 to indicate a level 0 element
2254  if (elem->level() == 0)
2255  {
2256  conn.push_back(-1);
2257  conn.push_back(-1);
2258  }
2259  else
2260  {
2261  conn.push_back(elem->parent()->id());
2262  conn.push_back(elem->parent()->which_child_am_i(elem));
2263  }
2264 #else
2265  conn.push_back (-1);
2266  conn.push_back (-1);
2267 #endif
2268 
2269  for (unsigned int n=0; n<elem->n_nodes(); n++)
2270  conn.push_back (elem->node(n));
2271 
2272  for (unsigned int n=0; n<elem->n_neighbors(); n++)
2273  {
2274  Elem *neigh = elem->neighbor(n);
2275  if (neigh)
2276  conn.push_back (neigh->id());
2277  else
2278  conn.push_back (-1);
2279  }
2280 
2281  elem->pack_indexing(std::back_inserter(conn));
2282 }
2283 
2284 
2285 
2287 {
2288 
2289  Elem *elem = Elem::build(this->type(),parent).release();
2290  libmesh_assert (elem);
2291 
2292 #ifdef LIBMESH_ENABLE_AMR
2293  if (this->level() != 0)
2294  {
2295  libmesh_assert(parent);
2296  parent->add_child(elem, this->which_child_am_i());
2297  libmesh_assert_equal_to (parent->type(), elem->type());
2298  libmesh_assert_equal_to (parent->child(this->which_child_am_i()), elem);
2299  }
2300 #endif
2301 
2302  // Assign the refinement flags and levels
2303 #ifdef LIBMESH_ENABLE_AMR
2304  elem->set_p_level(this->p_level());
2305  elem->set_refinement_flag(this->refinement_flag());
2307  libmesh_assert_equal_to (elem->level(), this->level());
2308 
2309  // If this element definitely should have children, assign
2310  // remote_elem for now; later unpacked elements may overwrite that.
2311  if (!elem->active())
2312  for (unsigned int c=0; c != elem->n_children(); ++c)
2313  elem->add_child(const_cast<RemoteElem*>(remote_elem), c);
2314 #endif
2315 
2316  // Assign the IDs
2317  elem->subdomain_id() = this->subdomain_id();
2318  elem->processor_id() = this->processor_id();
2319  elem->set_id() = this->id();
2320 
2321  // Assign the connectivity
2322  libmesh_assert_equal_to (elem->n_nodes(), this->n_nodes());
2323 
2324  for (unsigned int n=0; n<elem->n_nodes(); n++)
2325  elem->set_node(n) = mesh.node_ptr (this->node(n));
2326 
2327  // Assign the connectivity
2328  libmesh_assert_equal_to (elem->n_neighbors(), this->n_neighbors());
2329 
2330  for (unsigned int n=0; n<elem->n_neighbors(); n++)
2331  {
2332  dof_id_type neighbor_id = this->neighbor(n);
2333 
2334  // We should only be unpacking elements sent by their owners,
2335  // and their owners should know all their neighbors
2336  libmesh_assert_not_equal_to (neighbor_id, remote_elem->id());
2337 
2338  if (neighbor_id == DofObject::invalid_id)
2339  continue;
2340 
2341  Elem *neigh = mesh.query_elem(neighbor_id);
2342  if (!neigh)
2343  {
2344  elem->set_neighbor(n, const_cast<RemoteElem*>(remote_elem));
2345  continue;
2346  }
2347 
2348  // We never have neighbors more refined than us
2349  libmesh_assert_less_equal (neigh->level(), elem->level());
2350 
2351  // We never have subactive neighbors of non subactive elements
2352  libmesh_assert(!neigh->subactive() || elem->subactive());
2353 
2354  // If we have a neighbor less refined than us then it must not
2355  // have any more refined active descendants we could have
2356  // pointed to instead.
2357  libmesh_assert(neigh->level() == elem->level() ||
2358  neigh->active());
2359 
2360  elem->set_neighbor(n, neigh);
2361 
2362  // If neigh is at elem's level, then its family might have
2363  // remote_elem neighbor links which need to point to elem
2364  // instead, but if not, then we're done.
2365  if (neigh->level() != elem->level())
2366  continue;
2367 
2368  // What side of neigh is elem on? We can't use the usual Elem
2369  // method because we haven't finished restoring topology
2370  const AutoPtr<Elem> my_side = elem->side(n);
2371  unsigned int nn = 0;
2372  for (; nn != neigh->n_sides(); ++nn)
2373  {
2374  const AutoPtr<Elem> neigh_side = neigh->side(nn);
2375  if (*my_side == *neigh_side)
2376  break;
2377  }
2378 
2379  // elem had better be on *some* side of neigh
2380  libmesh_assert_less (nn, neigh->n_sides());
2381 
2382  // Find any elements that ought to point to elem
2383  std::vector<const Elem*> neigh_family;
2384 #ifdef LIBMESH_ENABLE_AMR
2385  if (!neigh->subactive())
2386  neigh->family_tree_by_side(neigh_family, nn);
2387 #else
2388  neigh_family.push_back(neigh);
2389 #endif
2390 
2391  // And point them to elem
2392  for (unsigned int i = 0; i != neigh_family.size(); ++i)
2393  {
2394  Elem* neigh_family_member = const_cast<Elem*>(neigh_family[i]);
2395 
2396  // The neighbor link ought to either be correct already or
2397  // ought to be to remote_elem
2398  libmesh_assert(neigh_family_member->neighbor(nn) == elem ||
2399  neigh_family_member->neighbor(nn) == remote_elem);
2400 
2401  neigh_family_member->set_neighbor(nn, elem);
2402  }
2403  }
2404 
2405  elem->unpack_indexing(this->indices());
2406 
2407  return elem;
2408 }
2409 
2410 
2411 
2412 unsigned int Elem::opposite_side(const unsigned int /*s*/) const
2413 {
2414  // If the subclass didn't rederive this, using it is an error
2415  libmesh_not_implemented();
2416 }
2417 
2418 
2419 
2420 unsigned int Elem::opposite_node(const unsigned int /*n*/,
2421  const unsigned int /*s*/) const
2422 {
2423  // If the subclass didn't rederive this, using it is an error
2424  libmesh_not_implemented();
2425 }
2426 
2427 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo