38 #include <boost/archive/binary_iarchive.hpp>
39 #include <boost/archive/binary_oarchive.hpp>
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>
53 #include <feel/feelcore/feel.hpp>
58 #include <feel/feelmesh/geoelement.hpp>
83 template<
typename Shape>
86 public VisitableBase<>,
88 public Elements<Shape>,
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 >
97 BOOST_STATIC_ASSERT( Shape::nDim == 3 );
102 static const uint16_type nDim = 3;
107 typedef typename VisitableBase<>::return_type return_type;
109 typedef VisitableBase<> super_visitable;
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;
120 typedef typename element_type::edge_permutation_type edge_permutation_type;
121 typedef typename element_type::face_permutation_type face_permutation_type;
123 typedef Points<3> super_points;
124 typedef typename super_points::points_type points_type;
125 typedef typename super_points::point_type point_type;
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;
132 typedef typename super_faces::face_iterator face_iterator;
133 typedef typename super_faces::face_const_iterator face_const_iterator;
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;
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;
146 typedef typename std::pair<size_type, size_type> edge_pair_type;
149 typedef boost::shared_ptr<self_type> self_ptrtype;
151 static const size_type SHAPE = Shape::Shape;
225 return super_elements::element_type::numLocalFaces;
233 return super_elements::element_type::numLocalEdges;
241 return super_elements::element_type::numLocalVertices;
249 return this->
faces().size();
257 return this->
edges().size();
266 return this->
points().size();
275 return M_e2e[e.id()][n];
299 virtual void setWorldComm( WorldComm
const& _worldComm )
301 this->setWorldCommMeshBase( _worldComm );
302 this->setWorldCommElements( _worldComm );
303 this->setWorldCommFaces( _worldComm );
304 this->setWorldCommEdges( _worldComm );
305 this->setWorldCommPoints( _worldComm );
312 virtual void clear();
315 FEELPP_DEFINE_VISITABLE();
329 FEELPP_ASSERT( 0 ).error(
"invalid call" );
362 friend class boost::serialization::access;
363 template<
class Archive>
364 void serialize( Archive & ar,
const unsigned int version )
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 );
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> );
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> );
398 boost::multi_array<element_edge_type,2> M_e2e;
401 template <
typename GEOSHAPE>
406 super_elements( worldComm ),
407 super_points( worldComm ),
408 super_faces( worldComm ),
409 super_edges( worldComm ),
413 template <
typename GEOSHAPE>
425 template <
typename GEOSHAPE>
429 template <
typename GEOSHAPE>
447 template <
typename GEOSHAPE>
453 this->
faces().clear();
454 this->
edges().clear();
456 M_e2e.resize( boost::extents[0][0] );
457 FEELPP_ASSERT(
isEmpty() ).error(
"all mesh containers should be empty after a clear." );
460 template <
typename GEOSHAPE>
463 std::vector<size_type>
const& cur, std::vector<uint32_type>& diff,
464 face_permutation_type& permutation, mpl::bool_<true> )
468 for ( uint16_type i = 0; i < def.size(); ++i )
469 diff[i] = def[i] - cur[2-i];
472 std::vector<uint32_type>::iterator _id_it = find( diff.begin(),
476 uint16_type pos = distance( diff.begin(), _id_it );
481 permutation = face_permutation_type::ROTATION_CLOCKWISE;
484 permutation = face_permutation_type::ROTATION_ANTICLOCK;
487 else if ( numZeros == 1 )
490 permutation = face_permutation_type::REVERSE_HYPOTENUSE;
493 permutation = face_permutation_type::REVERSE_HEIGHT;
496 permutation = face_permutation_type::REVERSE_BASE;
500 template <
typename GEOSHAPE>
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> )
506 std::vector<uint32_type>::iterator _id_it = find( diff.begin(),
510 uint16_type pos = distance( diff.begin(), _id_it );
515 permutation = face_permutation_type::SECOND_DIAGONAL;
518 permutation = face_permutation_type::PRINCIPAL_DIAGONAL;
521 else if ( numZeros == 0 )
523 if ( cur[0] == def[1] )
524 if ( cur[2] == def[3] )
525 permutation = face_permutation_type::REVERSE_BASE;
528 permutation = face_permutation_type::ROTATION_CLOCKWISE;
530 else if ( cur[0] == def[2] )
531 if ( cur[2] == def[0] )
532 permutation = face_permutation_type::ROTATION_ANTICLOCK;
535 permutation = face_permutation_type::REVERSE_HEIGHT;
538 permutation = face_permutation_type::ROTATION_TWICE_CLOCKWISE;
542 template <
typename GEOSHAPE>
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 );
552 for ( face_iterator elt_it = this->beginFace();
553 elt_it != this->endFace(); ++elt_it )
555 face_permutation_type permutation( face_permutation_type::IDENTITY );
561 for ( uint16_type i = 0; i < face_type::numVertices; ++i )
563 _left[i] = elt_it->element0().point( elt_it->element0().fToP( elt_it->pos_first(), i ) ).id();
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();
570 _diff[i] = _left[i] - _right[i];
573 uint16_type _numZeros = count( _diff.begin(), _diff.end(), uint32_type( 0 ) );
575 determineFacePermutation( _numZeros, _left, _right, _diff,
576 permutation, mpl::bool_<( SHAPE == SHAPE_TETRA )>() );
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(),
584 element_iterator iv, en;
587 for ( ; iv != en; ++iv )
589 for (
size_type j = 0; j < numLocalFaces(); j++ )
591 FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).warn(
"invalid element face check" );
595 DVLOG(2) <<
"[Mesh3D::updateFaces] element/face permutation : " << ti.elapsed() <<
"\n";
598 template <
typename GEOSHAPE>
603 std::map<std::set<int>,
size_type > _edges;
604 typename std::map<std::set<int>,
size_type >::iterator _edgeit;
606 M_e2e.resize( boost::extents[this->numElements()][this->numLocalEdges()] );
613 if ( ! this->
edges().empty() )
622 i1 = ( this->edge( j ).point( 0 ) ).
id();
623 i2 = ( this->edge( j ).point( 1 ) ).
id();
626 bool edgeinserted =
false;
628 boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
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" );
641 DVLOG(2) <<
"[Mesh3D::updateEdges] adding edges : " << ti.elapsed() <<
"\n";
645 edg.setProcessIdInPartition( this->worldComm().localRank() );
647 if ( this->
edges().empty() )
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 ) );
656 for ( ; ifa!=efa; ++ifa )
658 for ( uint16_type j = 0; j < face_type::numEdges; j++ )
661 i1 = bele.eToP( j, 0 );
662 i2 = bele.eToP( j, 1 );
664 i1 = ( ifa->point( i1 ) ).
id();
665 i2 = ( ifa->point( i2 ) ).
id();
668 bool edgeinserted =
false;
669 boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
680 edg.setId( _edgeit->second );
681 edg.setOnBoundary(
true, 0 );
683 for ( uint16_type k = 0; k < 2 + face_type::nbPtsPerEdge; k++ )
684 edg.setPoint( k, ifa->point( bele.eToP( j, k ) ) );
688 this->addEdge( edg );
691 if (!ifa->isGhostCell()) this->
edges().modify( this->edgeIterator( _edgeit->second ), Feel::detail::UpdateProcessId(ifa->processId()) );
696 DVLOG(2) <<
"[Mesh3D::updateEdges] adding edges : " << ti.elapsed() <<
"\n";
699 std::map<size_type, edge_pair_type> _oriented_edges;
700 typedef typename std::map<size_type, edge_pair_type>::iterator oe_iterator;
702 for ( element_iterator elt_it = this->beginElement();
703 elt_it != this->endElement(); ++elt_it )
707 for ( uint16_type j = 0; j < element_type::numEdges; ++j )
709 uint16_type j1 = ele.eToP( j, 0 );
710 uint16_type j2 = ele.eToP( j, 1 );
713 i1 = ( elt_it->point( j1 ) ).
id();
714 i2 = ( elt_it->point( j2 ) ).
id();
718 bool edgeinserted =
false;
719 boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
724 M_e2e[ vid ][ j] = boost::make_tuple( _edgeit->second, 1 );
731 FEELPP_ASSERT( _edgeit->second >= this->numEdges() )( _edgeit->second )( this->numEdges() ).error(
"invalid edge index" );
733 edg.setId( _edgeit->second );
736 edg.setOnBoundary(
false );
738 if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
739 edg.addElement( vid );
744 for ( uint16_type k = 0; k < 2 + element_type::nbPtsPerEdge; k++ )
745 edg.setPoint( k, elt_it->point( elt_it->eToP( j, k ) ) );
749 this->addEdge( edg );
754 if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
755 this->
edges().modify( this->edgeIterator( _edgeit->second ), [vid] ( edge_type& e ) { e.addElement( vid ); } );
759 if (!elt_it->isGhostCell()) this->
edges().modify( this->edgeIterator( _edgeit->second ), Feel::detail::UpdateProcessId(elt_it->processId()) );
762 detail::UpdateEdge<edge_type>( j, boost::cref( this->edge( _edgeit->second ) ) ) );
765 [j]( element_type
const& e ) { FEELPP_ASSERT( e.edgePtr( j ) )( e.id() )( j ).error(
"invalid edge in element" ); } );
768 for ( uint16_type j = 0; j < element_type::numFaces; ++j )
770 auto fit = this->
faces().iterator_to( elt_it->face(j));
771 for ( uint16_type e = 0; e < face_type::numEdges; ++e )
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); } );
780 DVLOG(2) <<
"[Mesh3D::updateEdges] updating element/edges : " << ti.elapsed() <<
"\n";
783 for ( element_iterator elt_it = this->beginElement();
784 elt_it != this->endElement(); ++elt_it )
792 i1 = ( elt_it->point( j1 ) ).
id();
793 i2 = ( elt_it->point( j2 ) ).
id();
797 bool edgeinserted =
false;
798 boost::tie( _edgeit, edgeinserted ) = _edges.insert( std::make_pair( s, next_edge ) );
803 edge_pair_type _current = std::make_pair( i1, i2 );
805 edge_permutation_type permutation( edge_permutation_type::IDENTITY );
807 oe_iterator _edge_it = _oriented_edges.find( _edgeit->second );
809 if ( _edge_it != _oriented_edges.end() )
811 edge_pair_type _default = _edge_it->second;
813 FEELPP_ASSERT( _default.first == _current.first ||
814 _default.first == _current.second ).error(
"invalid edge index" );
816 if ( _default.first != _current.first )
818 permutation = edge_permutation_type::REVERSE_PERMUTATION;
820 detail::UpdateEdgePermutation<edge_permutation_type>( j,
827 _oriented_edges.insert( std::make_pair( _edgeit->second, _current ) );
833 DVLOG(2) <<
"[Mesh3D::updateEdges] updating edges orientation : " << ti.elapsed() <<
"\n";
836 edge_iterator e_it = this->beginEdge();
837 edge_iterator e_en = this->endEdge();
839 for ( ; e_it!=e_en; ++e_it )
843 if ( e_it->numberOfElements() == 0 )
846 this->
edges().erase( e_it );
851 DVLOG(2) <<
"[Mesh3D::updateEdges] cleaning up edges : " << ti.elapsed() <<
"\n";
3D mesh class
Definition: mesh3d.hpp:84
elements_type const & elements() const
Definition: elements.hpp:355
element_edge_type const & localEdgeId(size_type const e, size_type const n) const
Definition: mesh3d.hpp:281
base mesh class
Definition: meshbase.hpp:67
~Mesh3D()
Definition: mesh3d.hpp:426
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
void updateEntitiesCoDimensionTwo()
Definition: mesh3d.hpp:600
size_type numLocalFaces() const
Definition: mesh3d.hpp:223
virtual void check() const =0
boost::tuple< size_type, size_type > face_processor_type
Definition: meshbase.hpp:83
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
Mesh3D(WorldComm const &worldComm=Environment::worldComm())
Definition: mesh3d.hpp:402
size_type numPoints() const
Definition: mesh3d.hpp:264
size_type numLocalEdges() const
Definition: mesh3d.hpp:231
size_type numLocalVertices() const
Definition: mesh3d.hpp:239
void renumber()
Definition: mesh3d.hpp:327
std::pair< element_iterator, element_iterator > elementsRange()
Definition: elements.hpp:452
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::pid_face_const_iterator, typename MeshTraits< MeshType >::pid_face_const_iterator > faces(MeshType const &mesh)
Definition: filters.hpp:933
boost::tuple< mpl::size_t< MESH_EDGES >, typename MeshTraits< MeshType >::pid_edge_const_iterator, typename MeshTraits< MeshType >::pid_edge_const_iterator > edges(MeshType const &mesh)
Definition: filters.hpp:1213
virtual bool isEmpty() const
Definition: elements.hpp:371
size_type numFaces() const
Definition: mesh3d.hpp:247
bool isEmpty() const
Definition: mesh3d.hpp:202
boost::tuple< size_type, int > element_edge_type
Definition: mesh3d.hpp:162
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Elements & operator=(Elements const &e)
Definition: elements.hpp:335
void updateEntitiesCoDimensionOnePermutation()
Definition: mesh3d.hpp:544
virtual void clear()
Definition: mesh3d.hpp:449
size_type numElements() const
Definition: mesh3d.hpp:215
element_edge_type const & localEdgeId(element_type const &e, size_type const n) const
Definition: mesh3d.hpp:272
size_type numEdges() const
Definition: mesh3d.hpp:255
WorldComm const & worldComm() const
Definition: meshbase.hpp:277