Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mesh3d.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-11-14
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007-2010 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __Mesh3D_H
31 #define __Mesh3D_H 1
32 
33 
34 #include <iomanip>
35 #include <fstream>
36 #include <cstdlib>
37 
38 #include <boost/archive/binary_iarchive.hpp>
39 #include <boost/archive/binary_oarchive.hpp>
40 
41 #include <boost/timer.hpp>
42 #include <boost/shared_ptr.hpp>
43 #include <boost/foreach.hpp>
44 #include <boost/multi_array.hpp>
45 #include <boost/multi_index_container.hpp>
46 #include <boost/multi_index/member.hpp>
47 #include <boost/multi_index/mem_fun.hpp>
48 #include <boost/multi_index/ordered_index.hpp>
49 #include <boost/numeric/ublas/io.hpp>
50 
51 
52 
53 #include <feel/feelcore/feel.hpp>
55 
57 
58 #include <feel/feelmesh/geoelement.hpp>
59 
61 #include <feel/feelmesh/faces.hpp>
62 #include <feel/feelmesh/edges.hpp>
63 #include <feel/feelmesh/points.hpp>
65 
66 namespace Feel
67 {
83 template<typename Shape>
84 class Mesh3D
85  :
86 public VisitableBase<>,
87 public MeshBase,
88 public Elements<Shape>,
89 public Points<3>,
90 public Faces<typename Shape::template shape<2>::type,
91  typename Elements<Shape>::element_type>,
92 public Edges<typename Shape::template shape<1>::type,
93  typename Faces<typename Shape::template shape<2>::type,
94  typename Elements<Shape>::element_type>::face_type >
95 {
96  // check at compilation time that the shape has indeed dimension 2
97  BOOST_STATIC_ASSERT( Shape::nDim == 3 );
98 
99  public:
100 
101 
102  static const uint16_type nDim = 3;
103 
107  typedef typename VisitableBase<>::return_type return_type;
108 
109  typedef VisitableBase<> super_visitable;
110  typedef MeshBase super;
111 
112  typedef Elements<Shape> super_elements;
113  typedef typename super_elements::elements_type elements_type;
114  typedef typename super_elements::element_type element_type;
115  typedef typename super_elements::element_iterator element_iterator;
116  typedef typename super_elements::element_const_iterator element_const_iterator;
117  typedef typename super_elements::update_element_neighbor_type update_element_neighbor_type;
118  typedef typename element_type::node_type node_type;
119 
120  typedef typename element_type::edge_permutation_type edge_permutation_type;
121  typedef typename element_type::face_permutation_type face_permutation_type;
122 
123  typedef Points<3> super_points;
124  typedef typename super_points::points_type points_type;
125  typedef typename super_points::point_type point_type;
126 
127  typedef Faces<typename Shape::template shape<2>::type,
128  typename super_elements::element_type> super_faces;
129  typedef typename super_faces::faces_type faces_type;
130  typedef typename super_faces::face_type face_type;
131 
132  typedef typename super_faces::face_iterator face_iterator;
133  typedef typename super_faces::face_const_iterator face_const_iterator;
134 
135  typedef typename super_faces::location_faces location_faces;
136  typedef typename super_faces::location_face_iterator location_face_iterator;
137  typedef typename super_faces::location_face_const_iterator location_face_const_iterator;
138 
139 
140  typedef Edges<typename Shape::template shape<1>::type,face_type> super_edges;
141  typedef typename super_edges::edges_type edges_type;
142  typedef typename super_edges::edge_type edge_type;
143  typedef typename super_edges::edge_iterator edge_iterator;
144  typedef typename super_edges::edge_const_iterator edge_const_iterator;
145 
146  typedef typename std::pair<size_type, size_type> edge_pair_type;
147 
148  typedef Mesh3D<Shape> self_type;
149  typedef boost::shared_ptr<self_type> self_ptrtype;
150 
151  static const size_type SHAPE = Shape::Shape;
152 
153  typedef typename super::face_processor_type face_processor_type;
154 
162  typedef boost::tuple<size_type, int> element_edge_type;
163 
165 
169 
173  Mesh3D( WorldComm const& worldComm = Environment::worldComm() );
174 
178  Mesh3D( Mesh3D const & m );
179 
183  ~Mesh3D();
184 
186 
190 
191  Mesh3D& operator=( Mesh3D const& m );
192 
194 
198 
202  bool isEmpty() const
203 {
204  return ( super_elements::isEmpty() &&
208  M_e2e.empty() );
209 }
210 
211 
216 {
217  return this->elements().size();
218 }
219 
224 {
225  return super_elements::element_type::numLocalFaces;
226 }
227 
232 {
233  return super_elements::element_type::numLocalEdges;
234 }
235 
240 {
241  return super_elements::element_type::numLocalVertices;
242 }
243 
248 {
249  return this->faces().size();
250 }
251 
256 {
257  return this->edges().size();
258 }
259 
260 
265 {
266  return this->points().size();
267 }
268 
272 element_edge_type const& localEdgeId( element_type const& e,
273  size_type const n ) const
274 {
275  return M_e2e[e.id()][n];
276 }
277 
282  size_type const n ) const
283 {
284  return M_e2e[e][n];
285 }
286 
288 
292 
293 
295 
299 virtual void setWorldComm( WorldComm const& _worldComm )
300 {
301  this->setWorldCommMeshBase( _worldComm );
302  this->setWorldCommElements( _worldComm );
303  this->setWorldCommFaces( _worldComm );
304  this->setWorldCommEdges( _worldComm );
305  this->setWorldCommPoints( _worldComm );
306 }
307 
312 virtual void clear();
313 
314 
315 FEELPP_DEFINE_VISITABLE();
316 
318 
319 
320 
321 protected:
322 
327 void renumber()
328 {
329  FEELPP_ASSERT( 0 ).error( "invalid call" );
330 }
331 
336 
341 
342 #if 0
343 
346 void updateFaces();
347 
351 void updateEdges();
352 
356 void check() const;
357 
358 #endif // 0
359 
360 private:
361 
362  friend class boost::serialization::access;
363  template<class Archive>
364  void serialize( Archive & ar, const unsigned int version )
365  {
366  ar & boost::serialization::base_object<super>( *this );
367  DVLOG(2) << "Serializing points\n";
368  ar & boost::serialization::base_object<super_points>( *this );
369  DVLOG(2) << "Serializing edges\n";
370  ar & boost::serialization::base_object<super_edges>( *this );
371  DVLOG(2) << "Serializing faces\n";
372  ar & boost::serialization::base_object<super_faces>( *this );
373  DVLOG(2) << "Serializing elements\n";
374  ar & boost::serialization::base_object<super_elements>( *this );
375  }
376 
377 
378 private:
379 
380 
384 void determineFacePermutation( uint16_type numZeros, std::vector<size_type> const& def,
385  std::vector<size_type> const& cur, std::vector<uint32_type>& diff,
386  face_permutation_type& permutation, mpl::bool_<true> );
387 
391 void determineFacePermutation( uint16_type numZeros, std::vector<size_type> const& def,
392  std::vector<size_type> const& cur, std::vector<uint32_type>& diff,
393  face_permutation_type& permutation, mpl::bool_<false> );
394 
398 boost::multi_array<element_edge_type,2> M_e2e;
399 };
400 
401 template <typename GEOSHAPE>
402 Mesh3D<GEOSHAPE>::Mesh3D( WorldComm const& worldComm )
403  :
404  super_visitable(),
405  super( worldComm ),
406  super_elements( worldComm ),
407  super_points( worldComm ),
408  super_faces( worldComm ),
409  super_edges( worldComm ),
410  M_e2e()
411 {}
412 
413 template <typename GEOSHAPE>
415  :
416  super_visitable(),
417  super( m ),
418  super_elements( m ),
419  super_points( m ),
420  super_faces( m ),
421  super_edges( m ),
422  M_e2e( m.M_e2e )
423 {}
424 
425 template <typename GEOSHAPE>
427 {}
428 
429 template <typename GEOSHAPE>
432 {
433  if ( this != &m )
434  {
435  super::operator=( m );
440 
441  M_e2e = m.M_e2e;
442  }
443 
444  return *this;
445 }
446 
447 template <typename GEOSHAPE>
448 void
450 {
451  this->elements().clear();
452  this->points().clear();
453  this->faces().clear();
454  this->edges().clear();
455 
456  M_e2e.resize( boost::extents[0][0] );
457  FEELPP_ASSERT( isEmpty() ).error( "all mesh containers should be empty after a clear." );
458 }
459 
460 template <typename GEOSHAPE>
461 void
462 Mesh3D<GEOSHAPE>::determineFacePermutation( uint16_type numZeros, std::vector<size_type> const& def,
463  std::vector<size_type> const& cur, std::vector<uint32_type>& diff,
464  face_permutation_type& permutation, mpl::bool_<true> )
465 {
466  if ( numZeros == 0 )
467  {
468  for ( uint16_type i = 0; i < def.size(); ++i )
469  diff[i] = def[i] - cur[2-i];
470  }
471 
472  std::vector<uint32_type>::iterator _id_it = find( diff.begin(),
473  diff.end(),
474  uint32_type( 0 ) );
475 
476  uint16_type pos = distance( diff.begin(), _id_it );
477 
478  if ( numZeros == 0 )
479  {
480  if ( pos == 0 )
481  permutation = face_permutation_type::ROTATION_CLOCKWISE;
482 
483  else
484  permutation = face_permutation_type::ROTATION_ANTICLOCK;
485  }
486 
487  else if ( numZeros == 1 )
488  {
489  if ( pos == 0 )
490  permutation = face_permutation_type::REVERSE_HYPOTENUSE;
491 
492  else if ( pos == 1 )
493  permutation = face_permutation_type::REVERSE_HEIGHT;
494 
495  else
496  permutation = face_permutation_type::REVERSE_BASE;
497  }
498 }
499 
500 template <typename GEOSHAPE>
501 void
502 Mesh3D<GEOSHAPE>::determineFacePermutation( uint16_type numZeros, std::vector<size_type> const& def,
503  std::vector<size_type> const& cur, std::vector<uint32_type>& diff,
504  face_permutation_type& permutation, mpl::bool_<false> )
505 {
506  std::vector<uint32_type>::iterator _id_it = find( diff.begin(),
507  diff.end(),
508  uint32_type( 0 ) );
509 
510  uint16_type pos = distance( diff.begin(), _id_it );
511 
512  if ( numZeros == 2 )
513  {
514  if ( pos == 0 )
515  permutation = face_permutation_type::SECOND_DIAGONAL;
516 
517  else
518  permutation = face_permutation_type::PRINCIPAL_DIAGONAL;
519  }
520 
521  else if ( numZeros == 0 )
522  {
523  if ( cur[0] == def[1] )
524  if ( cur[2] == def[3] )
525  permutation = face_permutation_type::REVERSE_BASE;
526 
527  else
528  permutation = face_permutation_type::ROTATION_CLOCKWISE;
529 
530  else if ( cur[0] == def[2] )
531  if ( cur[2] == def[0] )
532  permutation = face_permutation_type::ROTATION_ANTICLOCK;
533 
534  else
535  permutation = face_permutation_type::REVERSE_HEIGHT;
536 
537  else
538  permutation = face_permutation_type::ROTATION_TWICE_CLOCKWISE;
539  }
540 }
541 
542 template <typename GEOSHAPE>
543 void
545 {
546  boost::timer ti;
547  std::vector<size_type> _left( face_type::numVertices );
548  std::vector<size_type> _right( face_type::numVertices );
549  std::vector<uint32_type> _diff( face_type::numVertices );
550 
551  //determine permutation for the faces
552  for ( face_iterator elt_it = this->beginFace();
553  elt_it != this->endFace(); ++elt_it )
554  {
555  face_permutation_type permutation( face_permutation_type::IDENTITY );
556 
557  // if on boundary don't do anything
558  if ( elt_it->isOnBoundary() || ( elt_it->pos_second() == invalid_uint16_type_value ) )
559  continue;
560 
561  for ( uint16_type i = 0; i < face_type::numVertices; ++i )
562  {
563  _left[i] = elt_it->element0().point( elt_it->element0().fToP( elt_it->pos_first(), i ) ).id();
564 
565  uint16_type right_p = elt_it->element1().fToP( elt_it->pos_second(), i );
566  FEELPP_ASSERT( right_p >= 0 && right_p < elt_it->element1().numLocalPoints )( right_p )( elt_it->element1().numLocalPoints )
567  ( elt_it->pos_second() )( i ).error( "invalid point index" );
568  _right[i] = elt_it->element1().point( right_p ).id();
569 
570  _diff[i] = _left[i] - _right[i];
571  }
572 
573  uint16_type _numZeros = count( _diff.begin(), _diff.end(), uint32_type( 0 ) );
574 
575  determineFacePermutation( _numZeros, _left, _right, _diff,
576  permutation, mpl::bool_<( SHAPE == SHAPE_TETRA )>() );
577 
578  if ( permutation.value() != face_permutation_type::IDENTITY )
579  this->elements().modify( this->elementIterator( elt_it->ad_second(), elt_it->proc_second() ),
580  detail::UpdateFacePermutation<face_permutation_type>( elt_it->pos_second(),
581  permutation ) );
582  }
583 
584  element_iterator iv, en;
585  boost::tie( iv, en ) = this->elementsRange();
586 
587  for ( ; iv != en; ++iv )
588  {
589  for ( size_type j = 0; j < numLocalFaces(); j++ )
590  {
591  FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).warn( "invalid element face check" );
592  }
593  }
594 
595  DVLOG(2) << "[Mesh3D::updateFaces] element/face permutation : " << ti.elapsed() << "\n";
596 }
597 
598 template <typename GEOSHAPE>
599 void
601 {
602  boost::timer ti;
603  std::map<std::set<int>, size_type > _edges;
604  typename std::map<std::set<int>, size_type >::iterator _edgeit;
605  int next_edge = 0;
606  M_e2e.resize( boost::extents[this->numElements()][this->numLocalEdges()] );
607 
608  size_type vid, i1, i2;
609  element_type ele;
610  face_type bele;
611 
612  // First We check if we have already Edges stored
613  if ( ! this->edges().empty() )
614  {
615  // dump first the existing edges, to maintain the correct numbering
616  // if everything is correct, the numbering structure will reflect
617  // the actual edge numbering
618 
619  for ( size_type j = 0; j < this->edges().size(); ++j )
620  {
621  std::set<int> s;
622  i1 = ( this->edge( j ).point( 0 ) ).id();
623  i2 = ( this->edge( j ).point( 1 ) ).id();
624  s.insert( i1 );
625  s.insert( i2 );
626  bool edgeinserted = false;
627 
628  boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
629 
630  if ( edgeinserted )
631  ++next_edge;
632 
633  FEELPP_ASSERT( edgeinserted )( i1 )( i2 )(j)( this->edge( j ).id() )( _edgeit->second )( this->edge( j ).marker() )
634  ( this->edge( _edgeit->second ).point( 0 ).node() )( this->edge( _edgeit->second ).point( 1 ).node() )
635  ( this->edge( j ).point( 0 ).node() )( this->edge( j ).point( 1 ).node() ).error( "Two identical Edges stored in EdgeList" );
636  FEELPP_ASSERT( _edgeit->second == this->edge( j ).id() )( _edgeit->second )( this->edge( j ).id() ).error( "Edges in EdgeList have inconsistent id" );
637 
638  }
639  }
640 
641  DVLOG(2) << "[Mesh3D::updateEdges] adding edges : " << ti.elapsed() << "\n";
642  ti.restart();
643 
644  edge_type edg;
645  edg.setProcessIdInPartition( this->worldComm().localRank() );
646 
647  if ( this->edges().empty() )
648  {
649  // We want that the first edges be those on the boundary, in order to obey the paradigm for
650  // a Mesh3D
651  const location_faces& location_index=this->faces().template get<detail::by_location>();
652  location_face_iterator ifa;
653  location_face_iterator efa;
654  boost::tie( ifa, efa ) = location_index.equal_range( boost::make_tuple( true ) );
655 
656  for ( ; ifa!=efa; ++ifa )
657  {
658  for ( uint16_type j = 0; j < face_type::numEdges; j++ )
659  {
660  std::set<int> s;
661  i1 = bele.eToP( j, 0 );
662  i2 = bele.eToP( j, 1 );
663  // go to global
664  i1 = ( ifa->point( i1 ) ).id();
665  i2 = ( ifa->point( i2 ) ).id();
666  s.insert( i1 );
667  s.insert( i2 );
668  bool edgeinserted = false;
669  boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
670 
671  if ( edgeinserted )
672  ++next_edge;
673 
674  // reset the process id (edge not connected to an active elt take this value)
675  edg.setProcessId( invalid_uint16_type_value );
676 
677  if ( edgeinserted )
678  {
679  // set edge id
680  edg.setId( _edgeit->second );
681  edg.setOnBoundary( true, 0 );
682 
683  for ( uint16_type k = 0; k < 2 + face_type::nbPtsPerEdge; k++ )
684  edg.setPoint( k, ifa->point( bele.eToP( j, k ) ) );
685 
686  // TODO: should assocate a marker to the edge here ?
687  //edg.addElement( ifa->ad_first() );
688  this->addEdge( edg );
689  }
690  // set the process id from element (only active element)
691  if (!ifa->isGhostCell()) this->edges().modify( this->edgeIterator( _edgeit->second ), Feel::detail::UpdateProcessId(ifa->processId()) );
692  }
693  }
694  }
695 
696  DVLOG(2) << "[Mesh3D::updateEdges] adding edges : " << ti.elapsed() << "\n";
697  ti.restart();
698 
699  std::map<size_type, edge_pair_type> _oriented_edges;
700  typedef typename std::map<size_type, edge_pair_type>::iterator oe_iterator;
701 
702  for ( element_iterator elt_it = this->beginElement();
703  elt_it != this->endElement(); ++elt_it )
704  {
705  vid = elt_it->id();
706 
707  for ( uint16_type j = 0; j < element_type::numEdges; ++j )
708  {
709  uint16_type j1 = ele.eToP( j, 0 );
710  uint16_type j2 = ele.eToP( j, 1 );
711 
712  // go to global
713  i1 = ( elt_it->point( j1 ) ).id();
714  i2 = ( elt_it->point( j2 ) ).id();
715  std::set<int> s;
716  s.insert( i1 );
717  s.insert( i2 );
718  bool edgeinserted = false;
719  boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
720 
721  if ( edgeinserted )
722  ++next_edge;
723 
724  M_e2e[ vid ][ j] = boost::make_tuple( _edgeit->second, 1 );
725 
726  // reset the process id (edge not connected to an active elt take this value)
727  edg.setProcessId( invalid_uint16_type_value );
728 
729  if ( edgeinserted )
730  {
731  FEELPP_ASSERT( _edgeit->second >= this->numEdges() )( _edgeit->second )( this->numEdges() ).error( "invalid edge index" );
732  // set edge id
733  edg.setId( _edgeit->second );
734  // we have already inserted edges on the boundary so
735  // this one _is_ not on the boundary
736  edg.setOnBoundary( false );
737 
738  if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
739  edg.addElement( vid );
740 
741  // number of points on the edge is 2 (number of
742  // vertices) plus the number of points in the
743  // interior of the edge
744  for ( uint16_type k = 0; k < 2 + element_type::nbPtsPerEdge; k++ )
745  edg.setPoint( k, elt_it->point( elt_it->eToP( j, k ) ) );
746 
747  // TODO: should assocate a marker to the edge here ?
748  //edg.addElement( vid );
749  this->addEdge( edg );
750  }
751 
752  else
753  {
754  if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
755  this->edges().modify( this->edgeIterator( _edgeit->second ), [vid] ( edge_type& e ) { e.addElement( vid ); } );
756  }
757 
758  // set the process id from element (only active element)
759  if (!elt_it->isGhostCell()) this->edges().modify( this->edgeIterator( _edgeit->second ), Feel::detail::UpdateProcessId(elt_it->processId()) );
760 
761  this->elements().modify( elt_it,
762  detail::UpdateEdge<edge_type>( j, boost::cref( this->edge( _edgeit->second ) ) ) );
763 #if !defined(NDEBUG)
764  this->elements().modify( elt_it,
765  [j]( element_type const& e ) { FEELPP_ASSERT( e.edgePtr( j ) )( e.id() )( j ).error( "invalid edge in element" ); } );
766 #endif
767  }
768  for ( uint16_type j = 0; j < element_type::numFaces; ++j )
769  {
770  auto fit = this->faces().iterator_to( elt_it->face(j));
771  for ( uint16_type e = 0; e < face_type::numEdges; ++e )
772  {
773  auto const& elt_edge = elt_it->edge( elt_it->f2e( j, e ) );
774  this->faces().modify( fit,
775  [e,&elt_edge]( face_type& f ) { f.setEdge(e,elt_edge); } );
776  }
777  }
778  }
779 
780  DVLOG(2) << "[Mesh3D::updateEdges] updating element/edges : " << ti.elapsed() << "\n";
781  ti.restart();
782 
783  for ( element_iterator elt_it = this->beginElement();
784  elt_it != this->endElement(); ++elt_it )
785  {
786  for ( size_type j = 0; j < ( size_type )element_type::numEdges; ++j )
787  {
788  size_type j1 = ele.eToP( j, 0 );
789  size_type j2 = ele.eToP( j, 1 );
790 
791  // go to global
792  i1 = ( elt_it->point( j1 ) ).id();
793  i2 = ( elt_it->point( j2 ) ).id();
794  std::set<int> s;
795  s.insert( i1 );
796  s.insert( i2 );
797  bool edgeinserted = false;
798  boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
799 
800  if ( edgeinserted )
801  ++next_edge;
802 
803  edge_pair_type _current = std::make_pair( i1, i2 );
804 
805  edge_permutation_type permutation( edge_permutation_type::IDENTITY );
806 
807  oe_iterator _edge_it = _oriented_edges.find( _edgeit->second );
808 
809  if ( _edge_it != _oriented_edges.end() )
810  {
811  edge_pair_type _default = _edge_it->second;
812 
813  FEELPP_ASSERT( _default.first == _current.first ||
814  _default.first == _current.second ).error( "invalid edge index" );
815 
816  if ( _default.first != _current.first )
817  {
818  permutation = edge_permutation_type::REVERSE_PERMUTATION;
819  this->elements().modify( elt_it,
820  detail::UpdateEdgePermutation<edge_permutation_type>( j,
821  permutation ) );
822  }
823  }
824 
825  else
826  {
827  _oriented_edges.insert( std::make_pair( _edgeit->second, _current ) );
828  }
829 
830  }
831  }
832 
833  DVLOG(2) << "[Mesh3D::updateEdges] updating edges orientation : " << ti.elapsed() << "\n";
834  ti.restart();
835 #if 0
836  edge_iterator e_it = this->beginEdge();
837  edge_iterator e_en = this->endEdge();
838 
839  for ( ; e_it!=e_en; ++e_it )
840  {
841  // cleanup the edge data structure :
842 
843  if ( e_it->numberOfElements() == 0 )
844  {
845  // remove all edges that are not connected to any elements
846  this->edges().erase( e_it );
847  }
848 
849  }
850 
851  DVLOG(2) << "[Mesh3D::updateEdges] cleaning up edges : " << ti.elapsed() << "\n";
852 #endif
853  ti.restart();
854 }
855 
856 } // Feel
857 
858 #endif /* __Mesh3D_H */

Generated on Fri Oct 25 2013 14:24:19 for Feel++ by doxygen 1.8.4