29 #ifndef __MESHIMPL_HPP
30 #define __MESHIMPL_HPP 1
32 #include <boost/preprocessor/comparison/greater_equal.hpp>
44 template<
typename Shape,
typename T,
int Tag>
49 M_gm1( new gm1_type ),
53 M_tool_localization( new Localization() )
55 VLOG(2) <<
"[Mesh] constructor called\n";
57 template<
typename Shape,
typename T,
int Tag>
64 template<
typename Shape,
typename T,
int Tag>
68 flag_type theflag = -1;
69 if ( boost::any_cast<flag_type>( &__marker ) )
71 theflag = boost::any_cast<flag_type>( __marker);
73 else if ( boost::any_cast<int>( &__marker ) )
75 theflag = boost::any_cast<
int>( __marker);
77 else if ( boost::any_cast<size_type>( &__marker ) )
79 theflag = boost::any_cast<
size_type>( __marker);
81 else if ( boost::any_cast<uint16_type>( &__marker ) )
83 theflag = boost::any_cast<uint16_type>( __marker);
85 else if ( boost::any_cast<std::string>( &__marker ) )
87 theflag = this->markerName( boost::any_cast<std::string>( __marker) );
90 CHECK( theflag != -1 ) <<
"invalid flag type\n";
94 template<
typename Shape,
typename T,
int Tag>
98 VLOG(2) <<
"component MESH_RENUMBER: " << this->components().test( MESH_RENUMBER ) <<
"\n";
99 VLOG(2) <<
"component MESH_UPDATE_EDGES: " << this->components().test( MESH_UPDATE_EDGES ) <<
"\n";
100 VLOG(2) <<
"component MESH_UPDATE_FACES: " << this->components().test( MESH_UPDATE_FACES ) <<
"\n";
101 VLOG(2) <<
"component MESH_PARTITION: " << this->components().test( MESH_PARTITION ) <<
"\n";
103 if ( this->numElements() == 0 )
105 VLOG(2) <<
"No elements in Mesh? (with process rank " << this->worldComm().rank() <<
")\n";
106 if ( this->worldComm().localSize()==1 )
return;
111 VLOG(2) <<
"is already updated? : " << this->isUpdatedForUse() <<
"\n";
112 if ( !this->isUpdatedForUse() )
114 if ( this->components().test( MESH_RENUMBER ) )
118 VLOG(2) <<
"[Mesh::updateForUse] renumber : " << ti.elapsed() <<
"\n";
125 VLOG(2) <<
"Compute adjacency graph\n";
126 element_iterator iv, en;
127 std::map<size_type,boost::tuple<size_type, uint16_type, size_type> > f2e;
129 std::map<std::set<int>,
size_type > _faces;
130 typename std::map<std::set<int>,
size_type >::iterator _faceit;
135 for ( ; iv != en; ++iv )
137 for (
size_type j = 0; j < this->numLocalFaces(); j++ )
141 for (
int f = 0; f < face_type::numVertices; ++f )
144 s.insert( iv->point( j ).id() );
147 s.insert( iv->point( iv->fToP( j, f ) ).
id() );
150 bool faceinserted =
false;
151 boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
156 DVLOG(2) <<
"------------------------------------------------------------\n";
157 DVLOG(2) <<
"Element id: " << iv->id() <<
" local face id: " << j <<
" process id:" << iv->processId() <<
"\n";
161 DVLOG(2) <<
" new face " << _faceit->second <<
" is now in store with elt "
162 << f2e[_faceit->second].template get<0>() <<
" and local face id "
163 << f2e[_faceit->second].template get<1>() <<
"\n";
166 f2e[_faceit->second].template get<0>() = iv->id();
167 f2e[_faceit->second].template get<1>() = j;
173 DVLOG(2) <<
"old face " << _faceit->second <<
" was already in store with elt " << f2e[_faceit->second].template get<0>() <<
" and local face id " << f2e[_faceit->second].template get<1>() <<
"\n";
175 f2e[_faceit->second].template get<2>() = iv->id();
176 M_e2f[std::make_pair(iv->id(),j)]=boost::make_tuple( _faceit->second, f2e[_faceit->second].template get<0>() );
177 M_e2f[std::make_pair(f2e[_faceit->second].template get<0>(), f2e[_faceit->second].template get<1>())] =
178 boost::make_tuple( _faceit->second, iv->id() );
185 VLOG(2) <<
"Compute adjacency graph done in " << ti.elapsed() <<
"\n";
189 if ( this->components().test( MESH_PARTITION ) || ( this->worldComm().localSize() > 1 ) )
193 VLOG(2) <<
"[Mesh::updateForUse] partition time : " << ti1.elapsed() <<
"\n";
198 if ( this->components().test( MESH_UPDATE_FACES ) )
202 this->updateEntitiesCoDimensionOne();
204 this->updateEntitiesCoDimensionOnePermutation();
209 VLOG(2) <<
"[Mesh::updateForUse] update entities of codimension 1 : " << ti.elapsed() <<
"\n";
213 if ( this->components().test( MESH_UPDATE_EDGES ) )
218 this->updateEntitiesCoDimensionTwo();
219 VLOG(2) <<
"[Mesh::updateForUse] update edges : " << ti.elapsed() <<
"\n";
222 if ( this->components().test( MESH_UPDATE_FACES ) ||
223 this->components().test( MESH_UPDATE_EDGES )
228 this->setUpdatedForUse(
true );
232 element_iterator iv, en;
233 if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
236 for ( ; iv != en; ++iv )
238 this->
elements().modify( iv,
typename super_elements::ElementConnectPointToElement() );
243 auto pc = M_gm->preCompute( M_gm, M_gm->referenceConvex().vertices() );
244 auto pcf = M_gm->preComputeOnFaces( M_gm, M_gm->referenceConvex().barycenterFaces() );
248 for ( ; iv != en; ++iv )
251 [=,&pc,&pcf]( element_type& e )
253 for (
int i = 0; i < e.numPoints; ++i )
254 e.point( i ).addElement( e.id() );
255 e.setMeshAndGm(
this, M_gm, M_gm1 );
256 e.updateWithPc(pc, boost::ref( pcf) );
259 lambda::bind( &element_type::setMeshAndGm,
261 this, M_gm, M_gm1 ) );
264 lambda::bind( &element_type::updateWithPc,
265 lambda::_1, pc, boost::ref( pcf ) ) );
268 M_meas += iv->measure();
269 auto _faces = iv->faces();
274 for ( ; _faces.first != _faces.second; ++_faces.first )
275 if ( ( *_faces.first ) && ( *_faces.first )->isOnBoundary() )
276 M_measbdy += ( *_faces.first )->measure();
283 for ( ; iv != en; ++iv )
286 for(
auto _elt: iv->pointElementNeighborIds() )
289 meas += this->element( _elt ).measure();
291 this->
elements().modify( iv, [meas]( element_type& e ){ e.setMeasurePointElementNeighbors( meas ); } );
294 for (
auto itf = this->beginFace(), ite = this->endFace(); itf != ite; ++ itf )
296 this->
faces().modify( itf,[
this]( face_type& f ) { f.setMesh(
this ); } );
300 #if defined(FEELPP_ENABLE_MPI_MODE)
302 if ( this->components().test( MESH_UPDATE_FACES ) && this->worldComm().localSize()>1 )
304 this->updateEntitiesCoDimensionOneGhostCell();
308 if ( this->components().test( MESH_PROPAGATE_MARKERS ) )
309 propagateMarkers(mpl::int_<nDim>() );
315 M_gm->initCache(
this );
316 M_gm1->initCache(
this );
318 M_tool_localization->setMesh( this->shared_from_this(),
false );
322 VLOG(2) <<
"[Mesh::updateForUse] total time : " << ti.elapsed() <<
"\n";
324 template<
typename Shape,
typename T,
int Tag>
328 std::map<std::vector<double>,
size_type> mapDel;
329 for(
auto it = m.beginPoint(), en = m.endPoint();
334 if( it->isOnBoundary() )
336 for(
auto pit = this->beginPoint(), pen = this->endPoint();
340 if ( pit->isOnBoundary() )
342 if ( ublas::norm_2( it->node() - pit->node() ) < 1e-10 )
345 std::vector<double> ublasCopy(nDim);
346 std::copy(it->node().begin(), it->node().end(), ublasCopy.begin());
348 mapDel.insert(std::pair<std::vector<double>,
size_type>(ublasCopy, pit->id()));
355 this->addPoint( *it );
357 for(
auto it = m.beginFace(), en = m.endFace();it != en;++it )
360 if ( !it->isOnBoundary() )
367 for(
auto it = m.beginElement(), en = m.endElement();it != en;++it )
369 element_type e = *it;
370 for( uint16_type p = 0; p < e.numPoints; ++p ) {
371 if ( e.point(p).isOnBoundary() ) {
372 std::vector<double> ublasCopy(nDim);
373 std::copy(e.point(p).node().begin(), e.point(p).node().end(), ublasCopy.begin());
374 auto itMap = mapDel.find(ublasCopy);
375 if (itMap != mapDel.end() )
378 e.setPoint(p, this->point(itMap->second));
384 this->setUpdatedForUse(
false );
385 this->updateForUse();
388 template<
typename Shape,
typename T,
int Tag>
390 Mesh<Shape, T, Tag>::propagateMarkers( mpl::int_<2> )
393 std::for_each( this->beginFace(), this->endFace(),
394 [
this]( face_type
const& f )
396 if ( f.marker().isOff() )
400 for(
int i = 0; i < face_type::numPoints; ++i )
402 if ( f.point( i ).marker().isOff() )
405 this->
points().modify( this->
points().iterator_to( f.point(i) ),
406 [&f] ( point_type& p )
408 p.setMarker( f.marker().value() );
416 template<
typename Shape,
typename T,
int Tag>
418 Mesh<Shape, T, Tag>::propagateMarkers( mpl::int_<3> )
421 std::for_each( this->beginEdge(), this->endEdge(),
422 [
this]( edge_type
const& e )
424 if ( e.marker().isOff() )
427 for(
int i = 0; i < edge_type::numPoints; ++i )
429 if ( e.point( i ).marker().isOff() )
432 this->
points().modify( this->
points().iterator_to( e.point(i) ),
433 [&e] ( point_type& p )
435 p.setMarker( e.marker().value() );
442 std::for_each( this->beginFace(), this->endFace(),
443 [
this]( face_type
const& f )
445 if ( f.marker().isOff() )
449 for(
int i = 0; i < face_type::numPoints; ++i )
451 if ( f.point( i ).marker().isOff() )
454 this->
points().modify( this->
points().iterator_to( f.point(i) ),
455 [&f] ( point_type& p )
457 p.setMarker( f.marker().value() );
463 for(
int i = 0; i < face_type::numEdges; ++i )
465 if ( f.edge( i ).marker().isOff() )
468 this->
edges().modify( this->
edges().iterator_to( f.edge(i) ),
469 [&f] ( edge_type& e )
471 e.setMarker( f.marker().value() );
479 template<
typename Shape,
typename T,
int Tag>
489 typedef std::map<size_type,point_type> ptmap_type;
492 typedef typename ptmap_type::iterator ptmap_iterator;
493 typedef std::vector<size_type>::iterator nm_iterator;
496 for ( element_const_iterator elt = this->beginElement();
497 elt != this->endElement(); ++elt )
499 element_type
const& __element = *elt;
501 DVLOG(2) <<
"mesh::renumber] element id " << __element.id() <<
" proc " << __element.processId() <<
"\n";
503 for (
int i = 0; i < __element.nPoints(); ++i )
505 DVLOG(2) <<
"point id = " << __element.point( i ).id() <<
"\n";
510 for (
int i = 0; i < __element.nPoints(); ++i )
512 point_type __pt = __element.point( i );
520 __id = node_map[__true_id];
522 DVLOG_IF(2,(__id >= next_free_node )) <<
"next_free_node = " << next_free_node
523 <<
" point id = " << __id <<
"\n";
527 if ( __id == next_free_node )
529 node_map[ __true_id ] = next_free_node;
530 pt_map[next_free_node] = __pt;
531 pt_map[next_free_node].setId( next_free_node );
535 else if ( __id > next_free_node )
540 nm_iterator nm_it = std::find( node_map.begin(),
544 if ( nm_it != node_map.end() )
548 node_map[ next_free_node ] = __id;
552 node_map[ __true_id ] = next_free_node;
553 pt_map[next_free_node] = __pt;
554 pt_map[next_free_node].setId( next_free_node );
556 DVLOG(2) <<
"next_free_node = " << next_free_node
557 <<
" swapping point " << __true_id
558 <<
" and point " << next_free_node <<
"\n";
568 VLOG(2) <<
"[mesh::renumber] done collecting ids\n";
570 std::vector<size_type> check_id( node_map );
571 std::unique( check_id.begin(), check_id.end() );
573 FEELPP_ASSERT( check_id.size() == node_map.size() )( node_map.size() )( check_id.size() ).error(
"all ids must be unique" );
578 if ( Environment::numberOfProcessors() == 1 )
580 FEELPP_ASSERT( std::find( node_map.begin(),node_map.end(),
invalid_size_type_value ) == node_map.end() ).warn(
"invalid size_type value found as id " );
587 ptmap_iterator ptmapit = pt_map.begin();
588 ptmap_iterator ptmapen = pt_map.end();
590 for ( ; ptmapit != ptmapen; ++ptmapit )
592 point_iterator ptit = this->
points().find( ptmapit->first );
593 #if !defined( NDEBUG )
594 VLOG(2) <<
"[mesh::replace] replacing point " << ptit->id()
595 <<
" with " << ptmapit->second.id() <<
"\n";
597 bool __rep1 = this->
points().replace( ptit, ptmapit->second );
598 FEELPP_ASSERT( __rep1 )( __rep1 )( ptit->id() )( ptmapit->second.id() ) .warn(
"invalid point replacement" );
602 VLOG(2) <<
"[mesh::renumber] done replace point ids\n";
604 for ( element_iterator elt = this->beginElement();
605 elt != this->endElement(); ++elt )
607 element_type __element = *elt;
608 DVLOG(2) <<
"mesh::renumber] element id " << __element.id() <<
" proc " << __element.processId() <<
"\n";
611 for (
int i = 0; i < __element.nPoints(); ++i )
614 size_type __true_id =__element.point( i ).id();
615 auto& new_pt = this->point( node_map[__true_id] );
616 this->
elements().modify( elt,[i,&new_pt]( element_type& e ) { e.setPoint( i, new_pt ); } );
618 DVLOG(2) <<
"point id = " << __true_id <<
" node map " << node_map[__true_id] <<
" "
619 <<
"new point id = " << elt->point( i ).id() <<
"\n";
623 VLOG(2) <<
"[mesh::renumber] done replace point ids in elements\n";
625 for ( face_iterator elt = this->beginFace();
626 elt != this->endFace(); ++elt )
628 face_type __face = *elt;
629 DVLOG(2) <<
"face id: " << __face.id()
630 <<
" marker: " << __face.marker() <<
"\n";
633 for (
int i = 0; i < __face.nPoints(); ++i )
635 size_type __true_id =__face.point( i ).id();
636 this->
faces().modify( elt,
637 lambda::bind( &face_type::setPoint,
639 lambda::constant( i ),
640 boost::cref( this->point( node_map[__true_id] ) ) ) );
641 face_type __face2 = *elt;
643 DVLOG(2) <<
"id1= " << __face2.point( 0 ).id() <<
" id2= " << __face2.point( 1 ).id()<<
"\n";
644 DVLOG(2) <<
"point lid = " << i <<
" id = " << __true_id
645 <<
" nid = " << this->point( node_map[__true_id] ).id()
646 <<
" new point id = " << elt->point( i ).id() <<
"\n";
649 renumber( node_map, mpl::int_<nDim>() );
651 if ( Shape == SHAPE_TETRA && nOrder==1 )
656 checkLocalPermutation( mpl::bool_< ( Shape == SHAPE_TETRA ) >() );
665 template<
typename Shape,
typename T,
int Tag>
670 template<
typename Shape,
typename T,
int Tag>
676 template<
typename Shape,
typename T,
int Tag>
680 for (
auto elt = this->beginEdge();
681 elt != this->endEdge(); ++elt )
683 edge_type __edge = *elt;
685 DVLOG(2) <<
"edge id: " << __edge.id()
686 <<
" marker: " << __edge.marker() <<
"\n";
689 for (
int i = 0; i < __edge.nPoints(); ++i )
691 size_type __true_id =__edge.point( i ).id();
692 this->
edges().modify( elt,
693 lambda::bind( &edge_type::setPoint,
695 lambda::constant( i ),
696 boost::cref( this->point( node_map[__true_id] ) ) ) );
697 edge_type __edge2 = *elt;
699 DVLOG(2) <<
"renumber edge: id1= " << __edge2.point( 0 ).id() <<
" id2= " << __edge2.point( 1 ).id()<<
"\n";
700 DVLOG(2) <<
"renumber edge: point lid = " << i <<
" id = " << __true_id
701 <<
" nid = " << this->point( node_map[__true_id] ).id()
702 <<
" new point id = " << elt->point( i ).id() <<
"\n";
707 template<
typename Shape,
typename T,
int Tag>
711 for ( element_const_iterator elt = this->beginElement();
712 elt != this->endElement(); ++elt )
714 element_type __element = *elt;
716 int16_type n_swap = 0;
718 for ( uint16_type i=Shape::numVertices-1; i > 1; --i )
720 uint16_type min_index = 0;
722 for ( uint16_type j=1; j < __element.nPoints()+i-( Shape::numVertices-1 ); ++j )
724 if ( __element.point( j ).id() < __element.point( min_index ).id() )
728 if ( i != min_index )
730 __element.swapPoints( i, min_index );
731 n_swap += int16_type( 1 );
737 __element.swapPoints( 0,1 );
742 for ( uint16_type i=0; i < __element.nPoints(); ++i )
744 auto& new_pt = this->point( __element.point( i ).id() );
745 this->
elements().modify( elt, [i,&new_pt]( element_type& e ) { e.setPoint( i, new_pt ); } );
751 template<
typename Shape,
typename T,
int Tag>
756 updateEntitiesCoDimensionOne( mpl::bool_<true>() );
758 template<
typename Shape,
typename T,
int Tag>
762 template<
typename Shape,
typename T,
int Tag>
769 face.setProcessIdInPartition( this->worldComm().localRank() );
771 std::map<std::set<int>,
size_type > _faces;
772 typename std::map<std::set<int>,
size_type >::iterator _faceit;
777 if ( ! this->
faces().empty() )
787 face_iterator __it = this->beginFace();
788 face_iterator __en = this->endFace();
789 LOG(INFO) <<
"We have " << std::distance( __it, __en ) <<
" faces in the database";
790 auto itb = this->beginFaceOnBoundary();
791 auto enb = this->beginFaceOnBoundary();
792 LOG(INFO) <<
"We have " << std::distance( itb, enb ) <<
" faces on the boundary in the database";
793 for ( ; __it!=__en; )
797 for (
int f = 0; f < face_type::numVertices; ++f )
799 s.insert( __it->point( f ).id() );
802 bool faceinserted =
false;
803 boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
808 DVLOG_IF(2,faceinserted) <<
"added face with id " << __it->id () <<
"\n";
809 DVLOG_IF(2,!faceinserted) <<
"not added face with id " << __it->id ()
810 <<
" was already face with id = " << _faceit->second <<
"\n";
813 if ( faceinserted ==
false )
819 face_iterator __other = this->
faces().find( face_type( _faceit->second ) );
820 FEELPP_ASSERT( __other->id() != theid )
822 ( theid ).error(
"faces should have different ids " );
824 if ( __it->marker() != __other->marker() )
826 this->
faces().modify( __other, [__it]( face_type& f ) { f.setMarker2( __it->marker().value() ); } );
829 __it = this->eraseFace( __it );
839 face_type __f = *__it;
840 __f.setId( _faceit->second );
845 DVLOG(2) <<
"set face id " << __f.id()
846 <<
" iterator id = " << __it->id()
847 <<
" check id = " << _faceit->second <<
"\n";
849 this->
faces().replace( __it, __f );
855 VLOG(2) <<
"[Mesh::updateFaces] adding faces : " << ti.elapsed() <<
"\n";
858 VLOG(2) <<
"[Mesh::updateFaces] numLocalFaces : " << this->numLocalFaces() <<
"\n";
859 VLOG(2) <<
"[Mesh::updateFaces] face_type::numVertices : " << face_type::numVertices <<
"\n";
860 element_iterator iv, en;
863 for ( ; iv != en; ++iv )
865 element_type
const& __element = *iv;
869 for (
size_type j = 0; j < this->numLocalFaces(); j++ )
871 DVLOG(2) <<
"------------------------------------------------------------\n";
872 DVLOG(2) <<
"Element id: " << iv->id() <<
" local face id: " << j <<
"\n";
876 for (
int f = 0; f < face_type::numVertices; ++f )
878 uint16_type pt_localid = ( nDim==1 )?j:iv->fToP( j, f );
879 s.insert( iv->point( pt_localid ).id() );
880 VLOG(3) <<
"add point local id " << f <<
" to face " << j <<
" " << iv->fToP( j, f )
881 <<
" global id " << iv->point( pt_localid ).id() <<
"\n";
884 bool faceinserted =
false;
885 boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
890 DVLOG(2) <<
"------------------------------------------------------------\n";
891 DVLOG(2) <<
"Element id: " << iv->id() <<
" local face id: " << j <<
"\n";
895 DVLOG(2) <<
"creating the face:" << _faceit->second <<
"\n";
898 face.setId( _faceit->second );
900 if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
901 face.addElement( __element_id );
904 face.setProcessId( __element.processId() );
910 for (
size_type k = 0; k < face_type::numPoints; ++k )
911 face.setPoint( k, iv->point( ele.fToP( j, k ) ) );
914 face.setConnection0( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
915 face.setOnBoundary(
true, face_type::nDim );
918 bool inserted =
false;
920 boost::tie( __fit, inserted ) = this->addFace( face );
921 FEELPP_ASSERT( inserted && __fit != this->endFace() )
925 ( face.id() ).error(
"invalid face iterator" );
926 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
929 DVLOG(2) <<
"Adding [new] face info : \n";
930 DVLOG(2) <<
"element id: " << __element_id <<
"\n";
931 DVLOG(2) <<
"process id: " << __fit->processId() <<
"\n";
932 DVLOG(2) <<
"id: " << __fit->id() <<
"\n";
933 DVLOG(2) <<
"bdy: " << __fit->isOnBoundary() <<
"\n";
934 DVLOG(2) <<
"marker: " << __fit->marker() <<
"\n";
935 DVLOG(2) <<
"ad_first: " << __fit->ad_first() <<
"\n";
936 DVLOG(2) <<
"pos_first: " << __fit->pos_first() <<
"\n";
937 DVLOG(2) <<
"proc_first: " << __fit->proc_first() <<
"\n";
938 DVLOG(2) <<
"ad_second: " << __fit->ad_second() <<
"\n";
939 DVLOG(2) <<
"pos_second: " << __fit->pos_second() <<
"\n";
940 DVLOG(2) <<
"proc_second: " << __fit->proc_second() <<
"\n";
946 DVLOG(2) <<
"found the face:" << _faceit->second <<
" in element " << __element_id <<
" and local face: " << j <<
"\n";
949 face_iterator __fit = this->
faces().find( face_type( _faceit->second ) );
950 FEELPP_ASSERT( __fit != this->endFace() )( _faceit->second ).error(
"face is not in face container" );
953 face_type face = *__fit;
954 DVLOG(2) <<
"the face id :" << __fit->id() <<
"\n";
956 if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
957 face.addElement( __element_id );
960 if ( __fit->isConnectedTo0() && __fit->connection0().template get<0>() == 0 && ( __element.id() == __fit->ad_first() ) )
962 DVLOG(2) <<
"fixing connection 0 in face\n";
964 auto connect0 = __fit->connection0();
965 connect0.template get<0>() = boost::addressof( __element );
966 face.setConnection0( connect0 );
969 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
972 if ( __fit->isConnectedTo1() && __fit->connection1().template get<0>() == 0 && ( __element.id() == __fit->ad_second() ) )
974 DVLOG(2) <<
"fixing connection 1 in face\n";
977 auto connect1 = __fit->connection1();
978 connect1.template get<0>() = boost::addressof( __element );
979 face.setConnection1( connect1 );
982 element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
983 this->
elements().modify( elt1, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
984 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
987 if ( __fit->isConnectedTo0() && __fit->isConnectedTo1() )
989 DVLOG(2) <<
"internal face, fixing process id if necessary\n";
990 if ( !iv->facePtr( j ) )
991 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
992 FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).warn(
"invalid element face error" );
993 FEELPP_ASSERT( face.isConnectedTo0() && face.isConnectedTo1() )
994 ( face.isConnectedTo0() )( face.isConnectedTo1() ).error (
"inconsistent data structure" );
995 if ( face.processId()!=this->worldComm().localRank() )
997 if ( ( face.element0().processId()==this->worldComm().localRank() ) ||
998 ( face.element1().processId()==this->worldComm().localRank() ) )
999 face.setProcessId( this->worldComm().localRank() );
1007 if ( !__fit->isConnectedTo0() )
1009 DVLOG(2) <<
"[updateFaces][boundary] element: " << __element_id
1010 <<
" face: " << j <<
" id: " << _faceit->second <<
"\n";
1013 face.setConnection0( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
1015 face.setProcessId( __element.processId() );
1016 face.setOnBoundary(
true, face_type::nDim );
1020 this->
faces().replace( __fit, face );
1022 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1024 DVLOG(2) <<
"adding [!isConnectedTo0] face info : \n";
1025 DVLOG(2) <<
"id: " << __fit->id() <<
"\n";
1026 DVLOG(2) <<
"process id: " << __fit->processId() <<
"\n";
1027 DVLOG(2) <<
"bdy: " << __fit->isOnBoundary() <<
"\n";
1028 DVLOG(2) <<
"marker: " << __fit->marker() <<
"\n";
1029 DVLOG(2) <<
"ad_first: " << __fit->ad_first() <<
"\n";
1030 DVLOG(2) <<
"pos_first: " << __fit->pos_first() <<
"\n";
1031 DVLOG(2) <<
"proc_first: " << __fit->proc_first() <<
"\n";
1032 DVLOG(2) <<
"ad_second: " << __fit->ad_second() <<
"\n";
1033 DVLOG(2) <<
"pos_second: " << __fit->pos_second() <<
"\n";
1034 DVLOG(2) <<
"proc_second: " << __fit->proc_second() <<
"\n";
1035 DVLOG(2) <<
"element process id: " << iv->processId() <<
"\n";
1039 else if ( !__fit->isConnectedTo1() )
1044 this->
faces().modify( __fit,
1045 detail::UpdateFaceConnection1<typename face_type::element_connectivity_type>( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) ) );
1048 FEELPP_ASSERT( __fit->isConnectedTo0() && __fit->isConnectedTo1() )
1049 ( __fit->isConnectedTo0() )( __fit->isConnectedTo1() ).error(
"invalid face connection" );
1051 detail::UpdateFaceConnection1<typename face_type::element_connectivity_type> update1( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
1053 face.setOnBoundary(
false );
1056 if ( face.processId()!=this->worldComm().localRank() )
1058 if ( ( face.element0().processId()==this->worldComm().localRank() ) || ( face.element1().processId()==this->worldComm().localRank() ) )
1059 face.setProcessId( this->worldComm().localRank() );
1062 this->
faces().replace( __fit, face );
1065 element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
1066 this->
elements().modify( elt1, update_element_neighbor_type( __fit->pos_first(),
1067 __fit->ad_second() ) );
1068 this->
elements().modify( iv, update_element_neighbor_type( __fit->pos_second(),
1069 __fit->ad_first() ) );
1073 element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
1074 this->
elements().modify( elt1, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1075 this->
elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1078 DVLOG(2) <<
"adding face info : \n";
1079 DVLOG(2) <<
"id: " << __fit->id() <<
"\n";
1080 DVLOG(2) <<
"process id: " << __fit->processId() <<
"\n";
1081 DVLOG(2) <<
"bdy: " << __fit->isOnBoundary() <<
"\n";
1082 DVLOG(2) <<
"marker: " << __fit->marker() <<
"\n";
1083 DVLOG(2) <<
"ad_first: " << __fit->ad_first() <<
"\n";
1084 DVLOG(2) <<
"pos_first: " << __fit->pos_first() <<
"\n";
1085 DVLOG(2) <<
"proc_first: " << __fit->proc_first() <<
"\n";
1086 DVLOG(2) <<
"ad_second: " << __fit->ad_second() <<
"\n";
1087 DVLOG(2) <<
"pos_second: " << __fit->pos_second() <<
"\n";
1088 DVLOG(2) <<
"proc_second: " << __fit->proc_second() <<
"\n";
1089 DVLOG(2) <<
"element1 process id: " << elt1->processId() <<
"\n";
1090 DVLOG(2) <<
"element2 process id: " << iv->processId() <<
"\n";
1093 CHECK( (__fit->processId() == __fit->proc_first()) ||
1094 (__fit->processId() == __fit->proc_second()) )
1095 <<
"invalid process id " << __fit->processId() <<
" with element proc first = " << __fit->proc_first()
1096 <<
" and element proc second " << __fit->proc_second();
1099 LOG_IF( WARNING, !iv->facePtr( j ) ) <<
"invalid element " << iv->id() <<
" with local face id " << j;
1104 face_iterator f_it = this->beginFace();
1105 face_iterator f_en = this->endFace();
1107 for ( ; f_it!=f_en; )
1111 if ( !f_it->isConnectedTo0() )
1113 DLOG(INFO) <<
"removing face id : " << f_it->id()
1114 <<
" marker : " << f_it->marker();
1116 f_it = this->
faces().erase( f_it );
1126 LOG(INFO) <<
"We have now " <<
nelements(
boundaryfaces(
this)) <<
" faces on the boundary in the database";
1127 VLOG(2) <<
"element/face connectivity : " << ti.elapsed() <<
"\n";
1131 template<
typename Shape,
typename T,
int Tag>
1133 Mesh<Shape, T, Tag>::modifyEdgesOnBoundary( face_iterator& it , mpl::bool_<true> )
1136 for (
int f = 0; f < face_type::numEdges; ++f )
1138 if ( it->edge( f ).isOnBoundary() == false )
1140 auto eit = this->edgeIterator( it->edge(f).id() );
1141 this->
edges().modify( eit,
1143 {e.setOnBoundary(
true, 1 );} );
1150 template<
typename Shape,
typename T,
int Tag>
1152 Mesh<Shape, T, Tag>::modifyEdgesOnBoundary( face_iterator& f, mpl::bool_<false> )
1156 template<
typename Shape,
typename T,
int Tag>
1158 Mesh<Shape, T, Tag>::modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<false> )
1162 template<
typename Shape,
typename T,
int Tag>
1164 Mesh<Shape, T, Tag>::modifyElementOnBoundaryFromEdge( element_iterator& iv, mpl::bool_<true> )
1167 bool isOnBoundary =
false;
1168 for (
size_type j = 0; j < iv->nEdges(); j++ )
1170 isOnBoundary |= iv->edge( j ).isOnBoundary();
1175 this->
elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, 1 ); } );
1176 return isOnBoundary;
1178 return isOnBoundary;
1181 template<
typename Shape,
typename T,
int Tag>
1183 Mesh<Shape, T, Tag>::updateOnBoundary()
1187 LOG(INFO) <<
"update boundary points...";
1188 for(
auto it = this->beginFace(), en = this->endFace(); it != en; ++ it )
1190 if ( it->isOnBoundary() == true )
1192 modifyEdgesOnBoundary( it, mpl::bool_<(nDim >= 3)>() );
1194 for (
int f = 0; f < face_type::numPoints; ++f )
1196 if ( it->point( f ).isOnBoundary() == false )
1198 auto pit = this->pointIterator( it->point(f).id() );
1199 this->
points().modify( pit,
1201 {p.setOnBoundary(
true, 0 );} );
1208 LOG(INFO) <<
"update boundary elements...";
1210 for (
auto iv = this->beginElement(), en = this->endElement();
1213 bool isOnBoundary =
false;
1216 for (
size_type j = 0; j < iv->nTopologicalFaces(); j++ )
1218 isOnBoundary |= iv->face( j ).isOnBoundary();
1223 VLOG(3) <<
"checking " << iv->nTopologicalFaces() <<
" faces, isOnBoundary: " << isOnBoundary <<
" face_type::nDim: " << face_type::nDim;
1224 this->
elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, face_type::nDim ); } );
1228 bool e_modified = modifyElementOnBoundaryFromEdge( iv, mpl::bool_<(nDim>=3)>() );
1236 for (
size_type j = 0; j < iv->nPoints(); j++ )
1238 isOnBoundary |= iv->point( j ).isOnBoundary();
1243 VLOG(3) <<
"checking " << iv->nPoints() <<
" points, isOnBoundary: " << isOnBoundary;
1244 this->
elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, 0 ); } );
1248 <<
" elements sharing a point, a edge or a face with the boundary in the database";
1251 VLOG(3) <<
"boundary element : " << e.id()
1252 <<
" entity on boundary max dim " << e.boundaryEntityDimension()
1253 <<
" process id : " << e.processId();
1257 template<
typename Shape,
typename T,
int Tag>
1259 Mesh<Shape, T, Tag>::removeFacesFromBoundary( std::initializer_list<uint16_type> markers )
1261 std::for_each( markers.begin(), markers.end(),
1262 [=]( uint16_type marker )
1265 LOG(INFO) <<
"removing " <<
nelements(range) <<
" faces marked " << marker <<
" from boundary faces\n";
1267 for(
auto it = range.template get<1>(), en = range.template get<2>(); it != en; ++it )
1269 if ( it->isOnBoundary() )
1271 DVLOG(3) <<
"removing face " << it->id() <<
"\n";
1272 auto it2 = this->
faces().template project<0>( it );
1273 this->
faces().modify( it2, []( face_type & f ) { f.setOnBoundary(
false ); } );
1279 CHECK( it2->isOnBoundary() ==false ) <<
" face should not be on the boundary anymore\n";
1280 CHECK( it->isOnBoundary() ==false ) <<
" face should not be on the boundary anymore\n";
1287 #if defined(FEELPP_ENABLE_MPI_MODE)
1288 template<
typename Shape,
typename T,
int Tag>
1292 typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_type;
1294 VLOG(2) <<
"[Mesh::updateEntitiesCoDimensionOneGhostCell] start on god rank "<< this->worldComm().godRank() <<
"\n";
1296 std::vector<int> nbMsgToSend( this->worldComm().localSize() );
1297 std::fill( nbMsgToSend.begin(),nbMsgToSend.end(),0 );
1299 std::vector< std::map<int,int> > mapMsg( this->worldComm().localSize() );
1305 for ( ; iv != en; ++iv )
1307 element_type
const& __element = *iv;
1308 const int IdProcessOfGhost = __element.processId();
1309 const size_type idInPartition = __element.idInOthersPartitions( IdProcessOfGhost );
1311 this->
elements().modify( this->elementIterator(iv->id(),IdProcessOfGhost) ,
typename super_elements::ElementGhostConnectPointToElement() );
1313 this->
elements().modify( this->elementIterator(iv->id(),IdProcessOfGhost) ,
typename super_elements::ElementGhostConnectEdgeToElement() );
1316 this->worldComm().localComm().send( IdProcessOfGhost , nbMsgToSend[IdProcessOfGhost], idInPartition );
1318 std::cout<<
"I am the proc" << this->worldComm().localRank()<<
" , I send to proc " << IdProcessOfGhost
1319 <<
" with tag "<< nbMsgToSend[IdProcessOfGhost]
1320 <<
" idSend " << idInPartition
1321 <<
" it_ghost->G() " << __element.G()
1325 mapMsg[IdProcessOfGhost].insert( std::make_pair( nbMsgToSend[IdProcessOfGhost],__element.id() ) );
1327 ++nbMsgToSend[IdProcessOfGhost];
1332 std::vector<int> nbMsgToRecv;
1333 mpi::all_to_all( this->worldComm().localComm(),
1339 for (
int proc=0; proc<this->worldComm().localSize(); ++proc )
1341 for (
int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
1345 this->worldComm().localComm().recv( proc, cpt, idRecv );
1347 std::cout<<
"I am the proc" << this->worldComm().localRank()<<
" I receive to proc " << proc
1348 <<
" with tag "<< cpt <<
" idRecv " << idRecv
1349 <<
" it_ghost->G() " << this->element( idRecv ).G()
1352 auto const& theelt = this->element( idRecv );
1356 resultghost_type idFacesWithBary(this->numLocalFaces(),boost::make_tuple(0, std::vector<double>(nRealDim) ));
1357 for (
size_type j = 0; j < this->numLocalFaces(); j++ )
1359 auto const& theface = theelt.face( j );
1360 idFacesWithBary[j].template get<0>() = theface.id();
1362 auto const& theGj = theface.vertices();
1365 typename Localization::matrix_node_type v( theGj.size1(), 1 );
1366 ublas::scalar_vector<T> avg( theGj.size2(), T( 1 ) );
1367 T n_val = int( theGj.size2() );
1369 for (
size_type i = 0; i < theGj.size1(); ++i )
1370 v( i, 0 ) = ublas::inner_prod( ublas::row( theGj, i ), avg )/n_val;
1372 auto baryFace = ublas::column( v,0 );
1373 #else // doesn't work
1375 auto baryFace = ublas::column( glas::average( theface.G() ),0 );
1379 for ( uint16_type comp = 0; comp<nRealDim; ++comp )
1381 idFacesWithBary[j].template get<1>()[comp]=baryFace[comp];
1386 resultghost_type idPointsWithNode(element_type::numLocalVertices,boost::make_tuple(0, std::vector<double>(nRealDim) ));
1387 for (
size_type j = 0; j < element_type::numLocalVertices; j++ )
1389 auto const& thepoint = theelt.point( j );
1390 idPointsWithNode[j].template get<0>() = thepoint.id();
1391 for ( uint16_type comp = 0; comp<nRealDim; ++comp )
1393 idPointsWithNode[j].template get<1>()[comp]=thepoint(comp);
1397 std::vector<resultghost_type> theresponse(2);
1398 theresponse[0] = idPointsWithNode;
1399 theresponse[1] = idFacesWithBary;
1403 this->worldComm().localComm().send( proc, cpt, theresponse );
1409 for (
int proc=0; proc<this->worldComm().localSize(); ++proc )
1411 for (
int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1415 std::vector< boost::tuple<size_type, std::vector<double> > > idFacesWithBaryRecv(this->numLocalFaces());
1416 this->worldComm().localComm().recv( proc, cpt, idFacesWithBaryRecv );
1418 typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_face_type;
1419 typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_point_type;
1420 boost::tuple<resultghost_point_type,resultghost_face_type> requestRecv;
1421 this->worldComm().localComm().recv( proc, cpt, requestRecv );
1422 auto const& idPointsWithNodeRecv = requestRecv.template get<0>();
1423 auto const& idFacesWithBaryRecv = requestRecv.template get<1>();
1425 std::vector<resultghost_type> requestRecv;
1426 this->worldComm().localComm().recv( proc, cpt, requestRecv );
1427 auto const& idPointsWithNodeRecv = requestRecv[0];
1428 auto const& idFacesWithBaryRecv = requestRecv[1];
1432 std::cout<<
"I am the proc " << this->worldComm().localRank()<<
" I receive to proc " << proc
1433 <<
" with tag "<< cpt << std::endl;
1435 auto const& theelt = this->element( mapMsg[proc][cpt],proc );
1438 for (
size_type j = 0; j < this->numLocalFaces(); j++ )
1440 auto const& idFaceRecv = idFacesWithBaryRecv[j].template get<0>();
1441 auto const& baryFaceRecv = idFacesWithBaryRecv[j].template get<1>();
1446 for ( uint16_type j2 = 0; j2 < this->numLocalFaces() && !hasFind; j2++ )
1449 auto const& thefacej2 = theelt.face( j2 );
1451 auto const& theGj2 = thefacej2.vertices();
1454 typename Localization::matrix_node_type v( theGj2.size1(), 1 );
1455 ublas::scalar_vector<T> avg( theGj2.size2(), T( 1 ) );
1456 T n_val = int( theGj2.size2() );
1458 for (
size_type i = 0; i < theGj2.size1(); ++i )
1459 v( i, 0 ) = ublas::inner_prod( ublas::row( theGj2, i ), avg )/n_val;
1461 auto baryFace = ublas::column( v,0 );
1462 #else // doesn't compile (I don't now why)
1463 auto const& thefacej2 =theelt.face( j2 );
1464 auto baryFace = ublas::column( glas::average( thefacej2.G() ),0 );
1468 for (uint16_type d=0;d<nRealDim;++d)
1470 find2 = find2 && ( std::abs( baryFace[d]-baryFaceRecv[d] )<1e-9 );
1472 if (find2) { hasFind=
true;jBis=j2; }
1476 CHECK ( hasFind ) <<
"[mesh::updateEntitiesCoDimensionOneGhostCell] : invalid partitioning data, ghost face cells are not available\n";
1479 auto face_it = this->faceIterator( theelt.face( jBis ).id() );
1481 this->
faces().modify( face_it, detail::updateIdInOthersPartitions( proc, idFaceRecv ) );
1486 for (
size_type j = 0; j < element_type::numLocalVertices; j++ )
1488 auto const& idPointRecv = idPointsWithNodeRecv[j].template get<0>();
1489 auto const& nodePointRecv = idPointsWithNodeRecv[j].template get<1>();
1493 for ( uint16_type j2 = 0; j2 < element_type::numLocalVertices && !hasFind; j2++ )
1495 auto const& thepointj2 = theelt.point( j2 );
1498 for (uint16_type d=0;d<nRealDim;++d)
1500 find2 = find2 && ( std::abs( thepointj2(d)-nodePointRecv[d] )<1e-9 );
1502 if (find2) { hasFind=
true;jBis=j2; }
1505 CHECK ( hasFind ) <<
"[mesh::updateEntitiesCoDimensionOneGhostCell] : invalid partitioning data, ghost point cells are not available\n";
1507 auto point_it = this->pointIterator( theelt.point( jBis ).id() );
1509 this->
points().modify( point_it, detail::updateIdInOthersPartitions( proc, idPointRecv ) );
1524 #endif // defined(FEELPP_ENABLE_MPI_MODE)
1527 template<
typename Shape,
typename T,
int Tag>
1531 if ( nDim != nRealDim )
1533 #if !defined( NDEBUG )
1534 VLOG(2) <<
"[Mesh::check] numLocalFaces = " << this->numLocalFaces() <<
"\n";
1535 element_iterator iv = this->beginElementWithProcessId( this->worldComm().localRank() );
1536 element_iterator en = this->endElementWithProcessId( this->worldComm().localRank() );
1539 for ( ; iv != en; ++iv )
1541 element_type
const& __element = *iv;
1543 for (
size_type j = 0; j < this->numLocalFaces(); j++ )
1545 FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).error(
"invalid element face check" );
1546 VLOG(2) <<
"------------------------------------------------------------\n";
1547 VLOG(2) <<
"Element : " << iv->id() <<
" face lid: " << j <<
" face gid: " << iv->face( j ).id() <<
"\n";
1553 for ( uint16_type ms=0; ms < __element.nNeighbors(); ms++ )
1560 VLOG(2) <<
"[Mesh::check] element " << __element.id() <<
" number of neighbors: " << counter <<
"\n";
1561 FEELPP_ASSERT( counter >= 1 )( __element.id() )( __element.nNeighbors() )( counter ).warn(
"invalid neighboring data" );
1566 FEELPP_ASSERT( iv->edgePtr( j ) )( j )( iv->id() ).error(
"invalid element edge check" );
1567 VLOG(2) <<
"------------------------------------------------------------\n";
1568 VLOG(2) <<
"Element : " << iv->id() <<
" edge lid: " << j <<
" edge gid: " << iv->edge( j ).id() <<
"\n";
1577 typedef typename super::location_face_const_iterator location_face_const_iterator;
1578 location_face_const_iterator itf = this->beginFaceOnBoundary();
1579 location_face_const_iterator ite = this->endFaceOnBoundary();
1580 std::map<int,int> nf;
1581 for ( ; itf != ite; ++ itf )
1583 face_type
const& __face = *itf;
1584 DLOG_IF( WARNING, !__face.isConnectedTo0() && !__face.isConnectedTo1() ) <<
"face not connected to an element face:" << __face <<
"\n";
1585 if ( !__face.isConnectedTo0() && !__face.isConnectedTo1() )
1587 auto it = nf.find( (
int)__face.marker().value() );
1588 if ( it == nf.end() )
1589 nf[(
int)__face.marker().value()] = 1;
1593 DLOG_IF( WARNING, !__face.isConnectedTo0() && __face.isConnectedTo1() ) <<
"face connected to element 1 but not element 0. face:" << __face <<
"\n";
1594 DLOG_IF( WARNING, __face.isConnectedTo0() && __face.isConnectedTo1() )<<
"invalid boundary face (connected to 2 elements):" << __face <<
"\n";
1596 DLOG_IF( WARNING, __face.isConnectedTo0() && !__face.element( 0 ).facePtr( __face.pos_first() ) ) <<
"invalid face in element, face: " << __face <<
" in element " << __face.element(0) <<
"\n";
1599 auto itt = nf.begin();
1600 auto ent = nf.end();
1601 for(; itt != ent; ++itt )
1603 LOG(INFO ) <<
"face with marker " << itt->first <<
" not attached " << itt->second <<
"\n";
1607 template<
typename Shape,
typename T,
int Tag>
1613 if ( this->worldComm().localSize() == 1 )
1616 #ifdef FEELPP_HAS_MPI
1618 M_neighboring_processors.clear();
1621 Sphere bounding_sphere = processorBoundingSphere ( *
this, this->worldComm().localRank() );
1631 recv ( 4*this->worldComm().localSize(), 0 );
1633 send[0] = bounding_sphere.
center()( 0 );
1634 send[1] = bounding_sphere.
center()( 1 );
1635 send[2] = bounding_sphere.
center()( 2 );
1636 send[3] = bounding_sphere.
radius();
1638 MPI_Allgather ( &send[0], send.size(), MPI_DOUBLE,
1639 &recv[0], send.size(), MPI_DOUBLE,
1640 this->worldComm().localComm() );
1643 for (
unsigned int proc=0; proc<this->worldComm().localSize(); proc++ )
1645 const Point center ( recv[4*proc+0],
1649 const Real radius = recv[4*proc+3];
1651 const Sphere proc_sphere ( center, radius );
1653 if ( bounding_sphere.
intersects( proc_sphere ) )
1654 M_neighboring_processors.push_back( proc );
1658 VLOG(2) <<
"Processor " << this->worldComm().localRank() <<
" intersects:\n";
1660 for (
unsigned int p=0; p< M_neighboring_processors.size(); p++ )
1661 VLOG(2) <<
" - proc " << M_neighboring_processors[p] <<
"\n";
1667 template<
typename Shape,
typename T,
int Tag>
1671 bool mesh_well_oriented =
true;
1672 std::vector<size_type> list_of_bad_elts;
1674 for ( element_const_iterator elt = this->beginElement();
1675 elt != this->endElement(); ++elt )
1677 element_type
const& __element = *elt;
1679 if ( !__element.isAnticlockwiseOriented() )
1681 mesh_well_oriented =
false;
1682 list_of_bad_elts.push_back( elt->id() );
1686 if ( mesh_well_oriented )
1687 VLOG(2) <<
"Local numbering in the elements is OK . \n";
1691 std::for_each( list_of_bad_elts.begin(),
1692 list_of_bad_elts.end(),
1695 LOG_FIRST_N(WARNING,10) <<
"element is not anticlockwise oriented(wrong local numbering): " << e <<
"\n";
1702 template<
typename Shape,
typename T,
int Tag>
1706 element_type __element = *this->beginElement();
1708 if ( ! ( __element.isATriangleShape() ||
1709 __element.isATetrahedraShape() ) )
1712 static const uint16_type otn_triangle[ 10 ] =
1716 static const uint16_type otn_tetra[ 10 ] =
1718 1, 0, 2, 3, 4, 6, 5, 8, 7, 9
1721 for ( element_iterator elt = this->beginElement();
1722 elt != this->endElement(); ++elt )
1724 element_type __element = *elt;
1725 bool is_anticlockwise = __element.isAnticlockwiseOriented();
1728 FEELPP_ASSERT( is_anticlockwise ==
true )
1729 ( is_anticlockwise )
1731 ( __element.G() ).warn(
"invalid element permutation, will fix it" );
1735 if ( !is_anticlockwise )
1737 if ( __element.isATriangleShape() )
1738 __element.exchangePoints( otn_triangle );
1740 else if ( __element.isATetrahedraShape() )
1741 __element.exchangePoints( otn_tetra );
1747 ( __element.G() ).error(
"invalid element type" );
1748 throw std::logic_error(
"invalid element type" );
1751 is_anticlockwise = __element.isAnticlockwiseOriented();
1752 FEELPP_ASSERT( is_anticlockwise ==
true )
1753 ( is_anticlockwise )
1755 ( __element.G() ).error(
"invalid element permutation" );
1757 this->
elements().replace( elt, __element );
1763 template<
typename Shape,
typename T,
int Tag>
1768 VLOG(1) <<
"sending markername\n";
1770 VLOG(1) <<
"sending markername size: "<< M_markername.size() <<
"\n";
1771 BOOST_FOREACH(
auto m, M_markername )
1773 VLOG(1) <<
"sending key: "<< m.first <<
"\n";
1775 VLOG(1) <<
"sending value\n";
1780 template<
typename Shape,
typename T,
int Tag>
1784 VLOG(1) <<
"receiving markername\n";
1788 VLOG(1) <<
"receiving markername size: "<< s <<
"\n";
1789 for(
int i = 0; i < s; ++i )
1792 VLOG(1) <<
"receiving key\n";
1794 VLOG(1) <<
"receiving key:"<< k <<
"\n";
1796 VLOG(1) <<
"receiving value\n";
1798 VLOG(1) <<
"receiving value: "<< v[0] <<
","<< v[1] <<
"\n";
1807 template<
typename Shape,
typename T,
int Tag>
1808 typename Mesh<Shape, T, Tag>::element_iterator
1811 for(
int i = 0; i < element_type::numTopologicalFaces; ++i )
1813 auto fit = this->
faces().iterator_to( position->face(i) );
1814 auto modified = this->
faces().modify( fit, [position]( face_type& f ) { f.disconnect(*position); } );
1815 if ( fit->isOnBoundary() && !fit->isConnectedTo0() )
1817 std::cout <<
"erase boundary face...\n";
1818 this->eraseFace( fit );
1820 if ( !fit->isOnBoundary() && fit->isConnectedTo0() && !fit->isConnectedTo1() )
1822 std::cout <<
"found boundary face...\n";
1823 this->
faces().modify( fit, []( face_type& f ){ f.setOnBoundary(
true ); } );
1827 [] ( element_type& e )
1829 for (
int i = 0; i < e.numPoints; ++i )
1830 e.point( i ).elements().erase( e.id() );
1832 auto eit = this->
elements().erase( position );
1837 template<
typename Shape,
typename T,
int Tag>
1844 for(
auto pt_it = this->beginPoint(), pt_en = this->endPoint(); pt_it != pt_en; ++pt_it )
1846 std::vector<double> pts(3,0);
1847 pts[0] = pt_it->node()[0];
1848 if ( mesh_type::nRealDim >= 2 )
1849 pts[1] = pt_it->node()[1];
1850 if ( mesh_type::nRealDim >= 3 )
1851 pts[2] = pt_it->node()[2];
1852 M_enc_pts[pt_it->id()+1] = boost::make_tuple( pt_it->isOnBoundary(), pt_it->tags(), pts );
1859 auto face_it = allmarkedfaces.template get<1>();
1860 auto face_end = allmarkedfaces.template get<2>();
1864 GmshOrdering<element_type> ordering;
1866 GmshOrdering<face_type> ordering_face;
1869 for ( ; face_it != face_end; ++face_it )
1871 std::vector<int>
faces;
1872 faces.push_back( ordering_face.type() );
1873 faces.push_back( 4 + face_it->numberOfPartitions() );
1874 faces.push_back( face_it->marker().value() );
1875 faces.push_back( face_it->marker2().value() );
1876 faces.push_back( face_it->numberOfPartitions() );
1877 faces.push_back( face_it->processId() );
1878 for (
size_type i=0 ; i<face_it->numberOfNeighborPartitions(); ++i )
1879 faces.push_back( -( face_it->neighborPartitionIds()[i] ) );
1880 for ( uint16_type p=0; p<face_type::numPoints; ++p )
1881 faces.push_back( face_it->point( ordering_face.fromGmshId( p ) ).
id()+1 );
1883 M_enc_faces[face_it->id()] =
faces;
1887 auto eltOnProccess =
elements( *
this );
1888 auto elt_it = eltOnProccess.template get<1>();
1889 auto elt_en = eltOnProccess.template get<2>();
1891 for ( ; elt_it != elt_en; ++elt_it )
1893 std::vector<int> elts;
1894 elts.push_back( ordering.type() );
1895 elts.push_back( 4 + elt_it->numberOfPartitions() );
1896 elts.push_back( elt_it->marker().value() );
1897 elts.push_back( elt_it->marker2().value() );
1898 elts.push_back( elt_it->numberOfPartitions() );
1899 elts.push_back( elt_it->processId() );
1900 for (
size_type i=0 ; i<elt_it->numberOfNeighborPartitions(); ++i )
1901 elts.push_back( -( elt_it->neighborPartitionIds()[i] ) );
1902 for ( uint16_type p=0; p<element_type::numPoints; ++p )
1903 elts.push_back( elt_it->point( ordering.fromGmshId( p ) ).
id()+1 );
1905 M_enc_elts[elt_it->id()] = elts;
1910 template<
typename Shape,
typename T,
int Tag>
1915 std::vector<int> mapWorld(this->worldComm().size());
1916 for(
int cpu=0; cpu < this->worldComm().size(); ++cpu)
1917 mapWorld[cpu] = cpu;
1918 WorldComm worldcomm(mapWorld);
1925 this->setWorldComm(worldcomm.subWorldComm());
1927 LOG(INFO) <<
"decode= " << this->worldComm().size() <<
"\n" ;
1928 LOG(INFO) <<
"decode= " << this->worldComm().subWorldComm().localRank() <<
"\n";
1930 static const uint16_type npoints_per_face = ( face_type::numVertices*face_type::nbPtsPerVertex+
1931 face_type::numEdges*face_type::nbPtsPerEdge+
1932 face_type::numFaces*face_type::nbPtsPerFace );
1934 static const uint16_type npoints_per_element = element_type::numPoints;
1935 for(
auto pt_it = M_enc_pts.begin(), pt_en = M_enc_pts.end();
1936 pt_it != pt_en; ++pt_it )
1938 node_type __n( nRealDim );
1940 for ( uint16_type j = 0; j < nRealDim; ++j )
1941 __n[j] = pt_it->second.template get<2>()[j];
1943 point_type __pt( pt_it->first-1,__n, pt_it->second.template get<0>() );
1944 __pt.setOnBoundary( pt_it->second.template get<0>() );
1945 __pt.setTags( pt_it->second.template get<1>() );
1947 this->addPoint( __pt );
1949 for(
auto face_it = M_enc_faces.begin(), face_en = M_enc_faces.end();
1950 face_it != face_en; ++face_it )
1953 GmshOrdering<face_type> ordering;
1954 FEELPP_ASSERT( ordering.type() == face_it->second[0] )
1955 ( ordering.type() ) ( face_it->second[0] ).error(
"invalid mesh face type" );
1956 std::vector<int> tags( face_it->second[1] );
1957 for(
int i = 0; i < tags.size(); ++i ) tags[i] = face_it->second[2+i];
1959 pf.setId( this->numFaces() );
1960 pf.setProcessIdInPartition( this->worldComm().localRank() );
1961 pf.setProcessId( this->worldComm().localRank() );
1964 const int shift = face_it->second[1]+1;
1965 for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
1967 pf.setPoint( ordering.fromGmshId( jj ), this->point( face_it->second[shift+jj]-1 ) );
1969 this->addFace( pf );
1972 for(
auto elt_it = M_enc_elts.begin(), elt_en = M_enc_elts.end();
1973 elt_it != elt_en; ++elt_it )
1976 GmshOrdering<element_type> ordering;
1977 FEELPP_ASSERT( ordering.type() == elt_it->second[0] )
1978 ( ordering.type() ) ( elt_it->second[0] ).error(
"invalid mesh element type" );
1979 std::vector<int> tags( elt_it->second[1] );
1980 for(
int i = 0; i < tags.size(); ++i ) tags[i] = elt_it->second[2+i];
1982 pv.setProcessIdInPartition( this->worldComm().localRank() );
1983 pv.setProcessId( this->worldComm().localRank() );
1986 const int shift = elt_it->second[1]+1;
1987 for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1989 pv.setPoint( ordering.fromGmshId( jj ), this->point( elt_it->second[shift+jj]-1 ) );
1993 __idGmshToFeel=pv.id();
1994 auto theelt = mesh->elementIterator( pv.id(), pv.partitionId() );
1995 mesh->elements().modify( theelt, detail::updateIdInOthersPartitions( this->worldComm().localRank(), pv.id() ) );
1998 LOG(INFO) <<
"distance elts: "<< std::distance( this->beginElement(), this->endElement() ) <<
"\n";
1999 LOG(INFO) <<
"distance faces: "<< std::distance( this->beginFace(), this->endFace() ) <<
"\n";
2000 LOG(INFO) <<
"distance marker faces: "<< std::distance( this->beginFaceWithMarker(), this->endFaceWithMarker() ) <<
"\n";
2001 LOG(INFO) <<
"distance marker2 faces: "<< std::distance( this->beginFaceWithMarker2(), this->endFaceWithMarker2() ) <<
"\n";
2008 template<
typename Shape,
typename T,
int Tag1,
int Tag2=Tag1,
int TheTag=Tag1>
2009 boost::shared_ptr<Mesh<Shape, T, TheTag> >
2017 size_type shift_p = std::distance( m1->beginPoint(), m1->endPoint() );
2018 for(
auto it = m1->beginPoint(), en = m1->endPoint(); it != en; ++it )
2022 for(
auto it = m2->beginPoint(), en = m2->endPoint(); it != en; ++it )
2026 p.setId( shift_p+p.id() );
2029 typedef typename Mesh<Shape, T, TheTag>::face_iterator face_iterator;
2030 typedef typename Mesh<Shape, T, TheTag>::element_iterator element_iterator;
2031 typedef typename Mesh<Shape, T, TheTag>::face_type face_type;
2032 typedef typename Mesh<Shape, T, TheTag>::element_type element_type;
2034 auto addface = [&]( boost::shared_ptr<Mesh<Shape, T, TheTag> > m, face_iterator it,
size_type shift )
2038 pf.setOnBoundary(
true );
2039 const uint16_type npoints_per_face = ( face_type::numVertices*face_type::nbPtsPerVertex+
2040 face_type::numEdges*face_type::nbPtsPerEdge+
2041 face_type::numFaces*face_type::nbPtsPerFace );
2042 for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
2044 pf.setPoint( jj, m->point( shift+it->point(jj).id() ) );
2049 size_type shift_f = std::distance( m1->beginFace(), m1->endFace() );
2050 for(
auto it = m1->beginFace(), en = m1->endFace(); it != en; ++it )
2052 addface( m, it, 0 );
2054 for(
auto it = m2->beginFace(), en = m2->endFace(); it != en; ++it )
2057 addface( m, it, shift_p );
2060 auto addelement = [&]( boost::shared_ptr<Mesh<Shape, T, TheTag> > m, element_iterator it,
size_type shift_f,
size_type shift_p )
2062 element_type pf = *it;
2064 static const uint16_type npoints_per_element = element_type::numPoints;
2067 for ( uint16_type jj = 0; jj < pf.numLocalFaces; ++jj )
2069 pf.setFace( jj, m->face( shift_f+it->face(jj).id() ) );
2072 for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
2074 pf.setPoint( jj, m->point( shift_p+it->point(jj).id() ) );
2076 m->addElement( pf );
2078 for(
auto it = m1->beginElement(), en = m1->endElement(); it != en; ++it )
2080 addelement( m, it, 0, 0 );
2082 for(
auto it = m2->beginElement(), en = m2->endElement(); it != en; ++it )
2085 addelement( m, it, shift_f, shift_p );
2087 m->components().set ( MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
2092 template<
typename Shape,
typename T,
int Tag>
2094 Mesh<Shape, T, Tag>::Inverse::distribute(
bool extrapolation )
2096 typename self_type::element_iterator el_it;
2097 typename self_type::element_iterator el_en;
2098 boost::tie( boost::tuples::ignore, el_it, el_en ) =
Feel::elements( *M_mesh );
2100 typedef typename self_type::element_type element_type;
2101 typedef typename gm_type::template Context<vm::JACOBIAN|vm::KB|vm::POINT, element_type> gmc_type;
2102 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2105 typename gm_type::reference_convex_type refelem;
2106 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( M_mesh->gm(),
2107 refelem.points() ) );
2108 boost::unordered_map<size_type,bool> npt;
2111 M_ref_coords.clear();
2114 M_pts_cvx.resize( M_mesh->numElements() );
2116 KDTree::points_type boxpts;
2117 gmc_ptrtype __c(
new gmc_type( M_mesh->gm(),
2120 VLOG(2) <<
"[Mesh::Inverse] distribute mesh points ion kdtree\n";
2122 for ( ; el_it != el_en; ++el_it )
2125 __c->update( *el_it );
2126 gic_type gic( M_mesh->gm(), *el_it );
2130 bb.make( el_it->G() );
2132 for (
size_type k=0; k < bb.min.size(); ++k )
2138 DVLOG(2) <<
"G = " << el_it->G() <<
" min = " << bb.min <<
", max = " << bb.max <<
"\n";
2141 this->pointsInBox( boxpts, bb.min, bb.max );
2143 DVLOG(2) <<
"boxpts size = " << boxpts.size() <<
"\n";
2146 for (
size_type i = 0; i < boxpts.size(); ++i )
2148 size_type index = boost::get<1>( boxpts[i] );
2150 if ( ( !npt[index] ) || M_dist[index] < 0 )
2153 gic.setXReal( boost::get<0>( boxpts[i] ) );
2156 boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
2157 bool tobeadded = extrapolation || isin;
2159 DVLOG(2) <<
"i = " << i <<
" index = " << index <<
" isin = " << ( isin >= -1e-10 )
2160 <<
" xref = " << gic.xRef() <<
" xreal = " << boost::get<0>( boxpts[i] )
2161 <<
" tobeadded= " << tobeadded <<
" dist=" << dmin<<
"\n";
2164 if ( tobeadded && npt[index] )
2166 if ( dmin > M_dist[index] )
2167 M_pts_cvx[M_cvx_pts[index]].erase( index );
2175 M_ref_coords[index] = gic.xRef();
2176 M_dist[index] = dmin;
2177 M_cvx_pts[index]=el_it->id();
2178 M_pts_cvx[el_it->id()][index]=boost::get<2>( boxpts[i] );
2186 VLOG(2) <<
"[Mesh::Inverse] distribute mesh points in kdtree done\n";
2189 template<
typename Shape,
typename T,
int Tag>
2191 Mesh<Shape, T, Tag>::Localization::init()
2193 if ( !M_mesh )
return;
2195 DLOG_IF( WARNING, IsInit ==
false ) <<
"You have already initialized the tool of localization\n";
2199 M_geoGlob_Elts.clear();
2202 typename self_type::element_iterator el_it;
2203 typename self_type::element_iterator el_en;
2204 boost::tie( boost::tuples::ignore, el_it, el_en ) =
Feel::elements( *M_mesh );
2206 for ( ; el_it != el_en; ++el_it )
2208 for (
int i=0; i<el_it->nPoints(); ++i )
2210 if ( boost::get<1>( M_geoGlob_Elts[el_it->point( i ).id()] ).size()==0 )
2212 boost::get<0>( M_geoGlob_Elts[el_it->point( i ).id()] ) = el_it->point( i ).node();
2213 M_kd_tree->addPoint( el_it->point( i ).node(),el_it->point( i ).id() );
2216 boost::get<1>( M_geoGlob_Elts[el_it->point( i ).id()] ).push_back( el_it->id() );
2221 IsInitBoundaryFaces=
false;
2225 template<
typename Shape,
typename T,
int Tag>
2227 Mesh<Shape, T, Tag>::Localization::initBoundaryFaces()
2229 if ( !M_mesh )
return;
2231 #if !defined( NDEBUG )
2232 FEELPP_ASSERT( IsInitBoundaryFaces ==
false )
2233 ( IsInitBoundaryFaces ).warn(
"You have already initialized the tool of localization" );
2238 M_geoGlob_Elts.clear();
2241 typename self_type::location_face_iterator face_it;
2242 typename self_type::location_face_iterator face_en;
2245 for ( ; face_it != face_en; ++face_it )
2247 for (
int i=0; i<face_it->nPoints(); ++i )
2249 if ( face_it->isConnectedTo0() )
2251 if ( boost::get<1>( M_geoGlob_Elts[face_it->point( i ).id()] ).size()==0 )
2253 boost::get<0>( M_geoGlob_Elts[face_it->point( i ).id()] ) = face_it->point( i ).node();
2254 M_kd_tree->addPoint( face_it->point( i ).node(),face_it->point( i ).id() );
2256 boost::get<1>( M_geoGlob_Elts[face_it->point( i ).id()] ).push_back( face_it->element( 0 ).id() );
2261 IsInitBoundaryFaces=
true;
2268 template<
typename Shape,
typename T,
int Tag>
2269 boost::tuple<bool,typename Mesh<Shape, T, Tag>::node_type,
double>
2270 Mesh<Shape, T, Tag>::Localization::isIn(
size_type _id,
const node_type & _pt )
const
2277 auto const& elt = M_mesh->element( _id );
2279 if ( elt.isOnBoundary() )
2282 gmc_inverse_type gic( M_mesh->gm(), elt, this->mesh()->worldComm().subWorldCommSeq() );
2287 boost::tie( isin, dmin ) = M_refelem.isIn( gic.xRef() );
2292 gmc1_inverse_type gic( M_mesh->gm1(), elt, mpl::int_<1>(), this->mesh()->worldComm().subWorldCommSeq() );
2297 boost::tie( isin, dmin ) = M_refelem1.isIn( gic.xRef() );
2300 return boost::make_tuple(isin,x_ref,dmin);
2303 template<
typename Shape,
typename T,
int Tag>
2304 boost::tuple<uint16_type,std::vector<bool> >
2305 Mesh<Shape, T, Tag>::Localization::isIn( std::vector<size_type> _ids,
const node_type & _pt )
2307 typedef typename self_type::gm_type::reference_convex_type ref_convex_type;
2308 typedef typename self_type::gm1_type::reference_convex_type ref_convex1_type;
2310 uint16_type nbId = _ids.size();
2311 std::vector<bool> isin( _ids.size(),false );
2316 uint16_type nbIsIn=0;
2318 for ( uint16_type i = 0; i< nbId ; ++i )
2321 auto const& elt = M_mesh->element( _ids[i] );
2323 if ( elt.isOnBoundary() )
2326 gmc_inverse_type gic( M_mesh->gm(), elt );
2331 boost::tie( isin2, dmin ) = M_refelem.isIn( gic.xRef() );
2337 gmc1_inverse_type gic( M_mesh->gm1(), elt, mpl::int_<1>() );
2342 boost::tie( isin2, dmin ) = M_refelem1.isIn( gic.xRef() );
2345 if (isin[i]) ++nbIsIn;
2348 return boost::make_tuple( nbIsIn,isin );
2352 template<
typename Shape,
typename T,
int Tag>
2353 boost::tuple<bool, size_type, typename Mesh<Shape, T, Tag>::node_type>
2354 Mesh<Shape, T, Tag>::Localization::searchElement(
const node_type & p )
2357 #if !defined( NDEBUG )
2358 FEELPP_ASSERT( IsInit ==
true )
2359 ( IsInit ).warn(
"You don't have initialized the tool of localization" );
2361 bool isin=
false;
double dmin=0;
2363 size_type idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2365 std::list< std::pair<size_type, uint> > ListTri;
2366 this->searchInKdTree(p,ListTri);
2369 auto itLT=ListTri.begin();
2370 auto itLT_end=ListTri.end();
2372 #if !defined( NDEBUG )
2374 FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error(
" problem in list localization : is empty" );
2377 while ( itLT != itLT_end && !isin )
2383 boost::tie(isin,x_ref,dmin) = this->isIn(itLT->first,p);
2386 else idEltFound=itLT->first;
2392 if( this->doExtrapolation() )
2396 auto const& eltUsedForExtrapolation = this->mesh()->element(ListTri.begin()->first);
2397 gmc_inverse_type gic( this->mesh()->gm(), eltUsedForExtrapolation, this->mesh()->worldComm().subWorldCommSeq() );
2400 boost::tie(isin,idEltFound,x_ref) = boost::make_tuple(
true,eltUsedForExtrapolation.id(),gic.xRef());
2404 idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2410 return boost::make_tuple( isin, idEltFound, x_ref);
2414 template<
typename Shape,
typename T,
int Tag>
2415 boost::tuple<std::vector<bool>,
size_type>
2416 Mesh<Shape, T, Tag>::Localization::run_analysis(
const matrix_node_type & m,
2420 #if !defined( NDEBUG )
2421 FEELPP_ASSERT( IsInit ==
true )
2422 ( IsInit ).warn(
"You don't have initialized the tool of localization" );
2427 node_type x_ref;
double dmin=0;
2428 std::vector<bool> hasFindPts(m.size2(),
false);
2430 M_resultAnalysis.clear();
2432 bool doExtrapolationAtStart = this->doExtrapolation();
2433 auto nPtMaxNearNeighborAtStart = this->kdtree()->nPtMaxNearNeighbor();
2437 if ( doExtrapolationAtStart ) this->setExtrapolation(
false );
2440 auto currentEltHypothetical = eltHypothetical;
2441 for (
size_type i=0; i< m.size2(); ++i )
2443 bool testHypothetical_find =
false;
2447 boost::tie( testHypothetical_find,x_ref,dmin ) = this->isIn( currentEltHypothetical,ublas::column( m, i ) );
2449 if ( testHypothetical_find )
2451 cv_id = currentEltHypothetical;
2452 M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2458 boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2462 M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2463 currentEltHypothetical = cv_id;
2469 this->kdtree()->nbNearNeighbor(2*nPtMaxNearNeighborAtStart );
2471 boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2473 this->kdtree()->nbNearNeighbor(nPtMaxNearNeighborAtStart);
2477 M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2478 currentEltHypothetical = cv_id;
2481 else if (doExtrapolationAtStart)
2483 this->setExtrapolation(
true);
2485 boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2489 M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2490 currentEltHypothetical = cv_id;
2494 this->setExtrapolation(
false);
2506 this->setExtrapolation(doExtrapolationAtStart);
2508 return boost::make_tuple(hasFindPts,cv_id);
2514 template<
typename Shape,
typename T,
int Tag>
2515 boost::tuple<bool, std::list<boost::tuple<size_type, typename Mesh<Shape, T, Tag>::node_type> > >
2516 Mesh<Shape, T, Tag>::Localization::searchElements(
const node_type & p )
2519 #if !defined( NDEBUG )
2520 FEELPP_ASSERT( IsInit ==
true )
2521 ( IsInit ).warn(
"You don't have initialized the tool of localization" );
2526 std::list< std::pair<size_type, uint> > ListTri;
2527 searchInKdTree( p,ListTri );
2529 typename self_type::element_type elt;
2530 typename self_type::gm_type::reference_convex_type refelem;
2531 typename self_type::gm1_type::reference_convex_type refelem1;
2538 auto itLT=ListTri.begin();
2539 auto itLT_end=ListTri.end();
2541 #if !defined( NDEBUG )
2543 FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error(
" problem in list localization : is empty" );
2546 std::list<boost::tuple<size_type,node_type> > newlistelts;
2547 newlistelts.clear();
2549 bool finishSearch=
false;
2550 while ( itLT != itLT_end && !finishSearch )
2554 elt= M_mesh->element( itLT->first );
2556 if ( elt.isOnBoundary() )
2559 typename self_type::Inverse::gic_type gic( M_mesh->gm(), elt );
2566 boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
2572 typename self_type::Inverse::gic1_type gic( M_mesh->gm1(), elt,mpl::int_<1>() );
2579 boost::tie( isin, dmin ) = refelem1.isIn( gic.xRef() );
2583 boost::tie(isin,x_ref, dmin) = this->isIn(itLT->first,p);
2587 newlistelts.push_back( boost::make_tuple( itLT->first,x_ref ) );
2590 if (dmin>1e-7) finishSearch=
true;
2600 if ( !find ) std::cout <<
"\n WARNING EXTRAPOLATION IN SEARCHELEMENTS!!!"<<std::endl;
2603 return boost::make_tuple(
true,newlistelts );
2605 else if ( !find && !M_doExtrapolation )
2606 return boost::make_tuple(
false,newlistelts );
2611 itLT=ListTri.begin();
2612 elt= M_mesh->element( itLT->first );
2613 typename self_type::Inverse::gic_type gic( M_mesh->gm(), elt );
2619 newlistelts.push_back( boost::make_tuple( itLT->first,x_ref ) );
2621 return boost::make_tuple(
true,newlistelts );
2626 template<
typename Shape,
typename T,
int Tag>
2628 Mesh<Shape, T, Tag>::Localization::searchInKdTree(
const node_type & p,
2629 std::list< std::pair<size_type, uint> > & ListTri )
2632 M_kd_tree->search( p );
2635 typename KDTree::points_search_type ptsNN = M_kd_tree->pointsNearNeighbor();
2637 typename KDTree::points_search_const_iterator itNN = ptsNN.begin();
2638 typename KDTree::points_search_const_iterator itNN_end = ptsNN.end();
2640 #if !defined( NDEBUG )
2641 FEELPP_ASSERT( std::distance( itNN,itNN_end )>0 ).error(
"none Near Neighbor Points are find" );
2645 typename std::list<size_type>::iterator itL;
2646 typename std::list<size_type>::iterator itL_end;
2651 std::list< std::pair<size_type, uint> >::iterator itLT;
2652 std::list< std::pair<size_type, uint> >::iterator itLT_end;
2657 for ( ; itNN != itNN_end; ++itNN )
2659 itL= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).begin();
2660 itL_end= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).end();
2662 for ( ; itL != itL_end; ++itL )
2664 itLT=ListTri.begin();
2665 itLT_end=ListTri.end();
2668 while ( itLT != itLT_end && !find )
2670 if ( itLT->first == *itL ) find=
true;
2677 uint nb=itLT->second+1;
2679 ListTri.remove( *itLT );
2680 itLT=ListTri.begin();
2681 itLT_end=ListTri.end();
2684 while ( itLT != itLT_end && !find )
2686 if ( itLT->second < nb ) find=
true;
2691 ListTri.insert( itLT,std::make_pair( numEl,nb ) );
2694 else ListTri.push_back( std::make_pair( *itL,1 ) );
2700 template<
typename Shape,
typename T,
int Tag>
2701 boost::tuple<bool,typename Mesh<Shape, T, Tag>::node_type,
double>
2702 Mesh<Shape, T, Tag>::Localization::isIn(
size_type _id,
2703 const node_type & _pt,
2704 const matrix_node_type & setPoints,
2705 mpl::int_<1> )
const
2712 auto const& elt= M_mesh->element( _id );
2713 auto const& eltG = elt.G();
2716 std::vector<bool> find( setPoints.size2(),false );
2717 for (
size_type i=0; i< setPoints.size2() && isin; ++i )
2719 auto const& thePt = ublas::column( setPoints,i );
2721 for (
size_type j=0; j<eltG.size2(); ++j )
2723 auto ptjeltG = ublas::column( eltG,j );
2725 if ( ptjeltG.size()==1 )
2727 if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 )
2730 else if ( ptjeltG.size()==2 )
2732 if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 &&
2733 std::abs( thePt( 1 )-ptjeltG( 1 ) )<1e-8 )
2736 else if ( ptjeltG.size()==3 )
2738 if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 &&
2739 std::abs( thePt( 1 )-ptjeltG( 1 ) )<1e-8 &&
2740 std::abs( thePt( 2 )-ptjeltG( 2 ) )<1e-8 )
2752 boost::tie(isin2,x_ref,dmin) = this->isIn(_id,_pt);
2753 LOG_IF(WARNING, !isin2) <<
"Mesh::Localization::isIn<Conformal> : check fail -> maybe x_ref is not correct";
2756 return boost::make_tuple(isin,x_ref,dmin);
2761 template<
typename Shape,
typename T,
int Tag>
2762 boost::tuple<bool, size_type, typename Mesh<Shape, T, Tag>::node_type>
2763 Mesh<Shape, T, Tag>::Localization::searchElement(
const node_type & p,
2764 const matrix_node_type & setPoints,
2767 typename self_type::element_type elt;
2768 typename self_type::gm_type::reference_convex_type refelem;
2769 typename self_type::gm1_type::reference_convex_type refelem1;
2771 bool isin=
false,isin2=
false;
2774 size_type idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2776 std::list< std::pair<size_type, uint> > ListTri;
2777 searchInKdTree( p,ListTri );
2779 auto itLT=ListTri.begin();
2780 auto itLT_end=ListTri.end();
2782 #if !defined( NDEBUG )
2784 FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error(
" problem in list localization : is empty" );
2789 while ( itLT != itLT_end && !isin )
2791 boost::tie(isin,x_ref,dmin) = this->isIn(itLT->first,p,setPoints,mpl::int_<1>());
2794 if ( !isin ) ++itLT;
2795 else idEltFound=itLT->first;
2800 if( this->doExtrapolation() )
2802 std::cout <<
"WARNING EXTRAPOLATION for the point" << p << std::endl;
2804 auto const& eltUsedForExtrapolation = this->mesh()->element(ListTri.begin()->first);
2805 gmc_inverse_type gic( this->mesh()->gm(), eltUsedForExtrapolation, this->mesh()->worldComm().subWorldCommSeq() );
2808 boost::tie(isin,idEltFound,x_ref) = boost::make_tuple(
true,eltUsedForExtrapolation.id(),gic.xRef());
2812 idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2818 return boost::make_tuple( isin, idEltFound, x_ref);
2823 template<
typename Shape,
typename T,
int Tag>
2824 boost::tuple<std::vector<bool>,
size_type>
2825 Mesh<Shape, T, Tag>::Localization::run_analysis(
const matrix_node_type & m,
2827 const matrix_node_type & setPoints,
2830 #if !defined( NDEBUG )
2831 FEELPP_ASSERT( IsInit ==
true )
2832 ( IsInit ).warn(
"You don't have initialized the tool of localization" );
2839 std::vector<bool> hasFindPts(setPoints.size2(),
false);
2841 M_resultAnalysis.clear();
2842 auto currentEltHypothetical = eltHypothetical;
2843 for (
size_type i=0; i< m.size2(); ++i )
2846 bool testHypothetical_find =
false;
2849 boost::tie( testHypothetical_find,x_ref,dmin ) = this->isIn( currentEltHypothetical,ublas::column( m, i ),setPoints,mpl::int_<1>() );
2851 if ( testHypothetical_find )
2853 cv_id = currentEltHypothetical;
2854 M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2859 boost::tie( find_x, cv_id, x_ref ) = this->searchElement( ublas::column( m, i ),setPoints,mpl::int_<1>() );
2863 M_resultAnalysis[cv_id].push_back( boost::make_tuple( i,x_ref ) );
2865 currentEltHypothetical = cv_id;
2871 return boost::make_tuple(hasFindPts,cv_id);
2876 template<
typename Shape,
typename T,
int Tag>
2877 typename Mesh<Shape, T, Tag>::node_type
2878 Mesh<Shape, T, Tag>::Localization::barycenter()
const
2880 return this->barycenter(mpl::int_<nRealDim>());
2883 template<
typename Shape,
typename T,
int Tag>
2884 typename Mesh<Shape, T, Tag>::node_type
2885 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<1> )
const
2889 for (
size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2891 auto const& pt = this->kdtree()->points()[i].template get<0>();
2894 res(0)/=this->kdtree()->nPoints();
2897 template<
typename Shape,
typename T,
int Tag>
2898 typename Mesh<Shape, T, Tag>::node_type
2899 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<2> )
const
2903 for (
size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2905 auto const& pt = this->kdtree()->points()[i].template get<0>();
2906 res(0)+=pt(0);res(1)+=pt(1);
2908 res(0)/=this->kdtree()->nPoints();res(1)/=this->kdtree()->nPoints();
2911 template<
typename Shape,
typename T,
int Tag>
2912 typename Mesh<Shape, T, Tag>::node_type
2913 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<3> )
const
2916 res(0)=0;res(1)=0;res(2)=0;
2917 for (
size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2919 auto const& pt = this->kdtree()->points()[i].template get<0>();
2920 res(0)+=pt(0);res(1)+=pt(1);res(2)+=pt(2);
2922 res(0)/=this->kdtree()->nPoints();res(1)/=this->kdtree()->nPoints();res(2)/=this->kdtree()->nPoints();
2931 #if defined( FEELPP_INSTANTIATION_MODE )
2933 # define DIMS1 BOOST_PP_TUPLE_TO_LIST(1,(1))
2934 # define RDIMS1 BOOST_PP_TUPLE_TO_LIST(3,(1,2,3))
2935 # define ORDERS1 BOOST_PP_TUPLE_TO_LIST(5,(1,2,3,4,5))
2937 # define DIMS2 BOOST_PP_TUPLE_TO_LIST(1,(2))
2938 # define RDIMS2 BOOST_PP_TUPLE_TO_LIST(2,(2,3))
2939 # define ORDERS2 BOOST_PP_TUPLE_TO_LIST(5,(1,2,3,4,5))
2941 # define DIMS3 BOOST_PP_TUPLE_TO_LIST(1,(3))
2942 # define RDIMS3 BOOST_PP_TUPLE_TO_LIST(1,(3))
2943 # define ORDERS3 BOOST_PP_TUPLE_TO_LIST(4,(1,2,3,4))
2945 # define FACTORY_SIMPLEX(LDIM,LORDER,RDIM) template class Mesh<Simplex<LDIM, LORDER, RDIM> >;
2946 # define FACTORY_HYPERCUBE(LDIM,LORDER,RDIM) template class Mesh<Hypercube<LDIM, LORDER, RDIM> >;
2948 # define FACTORY_SIMPLEX_OP(_, GDO) FACTORY_SIMPLEX GDO
2949 # define FACTORY_HYPERCUBE_OP(_, GDO) FACTORY_HYPERCUBE GDO
2952 # define FACTORY_SIMPLEX_E(LDIM,LORDER,RDIM) extern template class Mesh<Simplex<LDIM, LORDER, RDIM> >;
2953 # define FACTORY_HYPERCUBE_E(LDIM,LORDER,RDIM) extern template class Mesh<Hypercube<LDIM, LORDER, RDIM> >;
2955 # define FACTORY_SIMPLEX_OP_E(_, GDO) FACTORY_HYPERCUBE GDO
2956 # define FACTORY_HYPERCUBE_OP_E(_, GDO) FACTORY_HYPERCUBE_E GDO
2958 #if !defined( FEELPP_MESH_IMPL_NOEXTERN )
2960 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS1, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS1 ), RDIMS1 ) )
2961 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS1, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS1 ), RDIMS1 ) )
2963 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS2, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS2 ), RDIMS2 ) )
2964 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS2, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS2 ), RDIMS2 ) )
2966 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS3, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS3 ), RDIMS3 ) )
2967 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS3, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS3 ), RDIMS3 ) )
2969 #endif // FEELPP_MESH_IMPL_NOEXTERN
2970 #endif // FEELPP_INSTANTIATION_MODE
2978 #endif // __MESHIMPL_HPP