31 #ifndef __ImporterGmsh_H
32 #define __ImporterGmsh_H 1
42 #include <boost/algorithm/string/trim.hpp>
49 void SwapBytes(
char *array,
int size,
int n);
65 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
69 x( Eigen::Vector3d::Zero() ),
70 uv( Eigen::Vector2d::Zero() ),
77 GMSHPoint(GMSHPoint
const& p )
83 parametric( p.parametric ),
99 is_on_processor( false ),
101 ghost_partition_id( -1 ),
114 std::vector<int>
const& _ghosts,
116 int _dom1,
int _dom2,
118 std::vector<int>
const& _indices,
119 int worldcommrank,
int worldcommsize)
125 numPartitions( _numPartitions ),
126 partition( _partition ),
128 is_on_processor( false ),
130 ghost_partition_id( partition ),
132 dom1(_dom1), dom2(_dom2),
133 numVertices( _numVertices ),
136 if ( worldcommsize == 1 )
138 is_on_processor =
true;
141 else if ( worldcommrank == partition )
143 is_on_processor =
true;
150 auto it = std::find( ghosts.begin(), ghosts.end(), worldcommrank );
151 if ( it != ghosts.end() )
153 is_on_processor =
true;
155 ghost_partition_id = partition;
160 GMSHElement( GMSHElement
const& g )
164 physical( g.physical ),
165 elementary( g.elementary ),
166 numPartitions( g.numPartitions ),
167 partition( g.partition ),
169 is_on_processor( g.is_on_processor ),
170 is_ghost( g.is_ghost ),
171 ghost_partition_id( g.ghost_partition_id ),
173 dom1( g.dom1 ), dom2( g.dom2 ),
174 numVertices( g.numVertices ),
178 bool isOnProcessor()
const {
return is_on_processor; }
179 bool isGhost()
const {
return is_ghost; }
180 int ghostPartitionId()
const {
return ghost_partition_id; }
182 template<
typename IteratorBegin,
typename IteratorEnd>
183 bool isIgnored(IteratorBegin it, IteratorEnd en )
const {
return std::find( it, en, physical ) != en; }
193 std::vector<int> ghosts;
194 bool is_on_processor;
196 int ghost_partition_id;
203 std::vector<int> indices;
226 template<
typename MeshType>
244 M_use_elementary_region_as_physical_region=param;
251 typedef typename super::mesh_type mesh_type;
252 typedef typename super::point_type point_type;
253 typedef typename super::node_type node_type;
254 typedef typename super::edge_type edge_type;
255 typedef typename super::face_type face_type;
256 typedef typename super::element_type element_type;
258 typedef typename mesh_type::face_iterator face_iterator;
264 static const uint16_type npoints_per_edge = ( edge_type::numVertices*edge_type::nbPtsPerVertex+
265 edge_type::numEdges*edge_type::nbPtsPerEdge+
266 edge_type::numFaces*edge_type::nbPtsPerFace );
268 static const uint16_type npoints_per_face = ( face_type::numVertices*face_type::nbPtsPerVertex+
269 face_type::numEdges*face_type::nbPtsPerEdge+
270 face_type::numFaces*face_type::nbPtsPerFace );
272 static const uint16_type npoints_per_element = element_type::numPoints;
279 ImporterGmsh( WorldComm
const& _worldcomm = Environment::worldComm() )
281 super( GMSH, _worldcomm ),
282 M_version( FEELPP_GMSH_FORMAT_VERSION ),
283 M_use_elementary_region_as_physical_region( false )
285 this->setIgnorePhysicalName(
"FEELPP_GMSH_PHYSICALNAME_IGNORED" );
289 explicit ImporterGmsh( std::string
const& _fname, std::string _version = FEELPP_GMSH_FORMAT_VERSION,
290 WorldComm
const& _worldcomm = Environment::worldComm() )
292 super( _fname, GMSH, _worldcomm ),
293 M_version( _version ),
294 M_use_elementary_region_as_physical_region( false )
296 this->setIgnorePhysicalName(
"FEELPP_GMSH_PHYSICALNAME_IGNORED" );
299 ImporterGmsh( ImporterGmsh
const & i )
302 M_version( i.M_version ),
303 M_use_elementary_region_as_physical_region( false ),
304 M_ignorePhysicalGroup( i.M_ignorePhysicalGroup ),
305 M_ignorePhysicalName( i.M_ignorePhysicalName )
307 this->setIgnorePhysicalName(
"FEELPP_GMSH_PHYSICALNAME_IGNORED" );
338 boost::tuple<bool,boost::tuple<bool,int> >
isElementOnProcessor( std::vector<int>
const& tags )
const;
346 void setVersion( std::string
const&
version )
351 void setIgnorePhysicalGroup(
int i )
353 M_ignorePhysicalGroup.insert( i );
355 void setIgnorePhysicalName( std::string s )
357 M_ignorePhysicalName.insert( s );
366 void visit( mesh_type* mesh );
378 void addPoint( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel );
379 void addPoint( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<1> );
380 void addPoint( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<2> );
381 void addPoint( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<3> );
383 void addEdge( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel );
384 void addEdge( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<1> );
385 void addEdge( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<2> );
386 void addEdge( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<3> );
388 void addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel );
389 void addFace( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<1> );
390 void addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<2> );
391 void addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<3> );
393 void addVolume( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel );
394 void addVolume( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<1> );
395 void addVolume( mesh_type* , Feel::detail::GMSHElement
const& ,
int & , mpl::int_<2> );
396 void addVolume( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & , mpl::int_<3> );
398 void updateGhostCellInfo( mesh_type* mesh, std::map<int,int>
const& __idGmshToFeel, std::map<
int,boost::tuple<int,int> >
const& __mapGhostElt );
403 std::string M_version;
404 std::map<int,int> M_n_vertices;
407 std::set<int> M_ignorePhysicalGroup;
408 std::set<std::string> M_ignorePhysicalName;
409 bool M_use_elementary_region_as_physical_region;
418 template<
typename MeshType>
420 ImporterGmsh<MeshType>::showMe()
const
422 DVLOG(2) <<
"[ImporterGmsh::showMe] npoints_per_element = " << npoints_per_element <<
"\n";
423 DVLOG(2) <<
"[ImporterGmsh::showMe] npoints_per_face = " << npoints_per_face <<
"\n";
424 DVLOG(2) <<
"[ImporterGmsh::showMe] npoints_per_edge = " << npoints_per_edge <<
"\n";
426 template<
typename MeshType>
427 boost::tuple<bool,boost::tuple<bool,int> >
430 bool is_element_on_processor =
false;
434 if ( tag.size() == 2 )
435 return boost::make_tuple(
true,
false );
439 int rank = this->worldComm().localRank();
440 auto it = std::find_if( tag.begin()+3, tag.end(), [&rank] (
int i )
445 if ( it != tag.end() )
446 is_element_on_processor =
true;
448 bool isGhostCell=
false;
454 idOfGhostCell=tag[3];
457 return boost::make_tuple( is_element_on_processor, boost::make_tuple( isGhostCell,idOfGhostCell ) );
459 template<
typename MeshType>
464 if ( this->version() !=
"1.0" &&
465 this->version() !=
"2.0" &&
466 this->version() !=
"2.1" &&
467 this->version() !=
"2.2" &&
468 this->version() != FEELPP_GMSH_FORMAT_VERSION )
469 throw std::logic_error(
"invalid gmsh file format version" );
471 DVLOG(2) <<
"visit(" << mesh_type::nDim <<
"D ) starts\n";
472 DVLOG(2) <<
"visit(" << mesh_type::nDim <<
"D ) filename = " << this->filename() <<
"\n";
474 std::ifstream __is ( this->filename().c_str() );
476 if ( !__is.is_open() )
478 std::ostringstream ostr;
479 LOG(ERROR) <<
"Invalid file name " << this->filename() <<
" (file not found)";
480 ostr <<
"Invalid file name " << this->filename() <<
" (file not found)\n";
481 throw std::invalid_argument( ostr.str() );
487 std::string theversion;
488 double version = 2.2;
490 bool binary =
false, swap =
false;
492 if ( ( ( this->version() ==
"2.0" ) ||
493 ( this->version() ==
"2.1" ) ||
494 ( this->version() ==
"2.2" ) ||
495 ( this->version() == FEELPP_GMSH_FORMAT_VERSION ) ) &&
496 std::string( __buf ) ==
"$MeshFormat" )
501 __is >> theversion >> format >> size;
502 LOG(INFO) <<
"GMSH mesh file version : " << theversion <<
" format: " << (format?
"binary":
"ascii") <<
" size of double: " << size <<
"\n";
503 CHECK( boost::lexical_cast<double>( theversion ) >= 2 )
504 <<
"Feel++ supports only Gmsh version >= 2.\n";
505 version = boost::lexical_cast<
double>( theversion );
509 CHECK( (c=__is.get()) ==
'\n' ) <<
"Invalid character " << c <<
" should be newline\n";
511 LOG(INFO) <<
"GMSH mesh is in binary format\n";
513 __is.read( (
char*)&one,
sizeof(
int) );
517 LOG(INFO) <<
"one before swap : " << one <<
"\n";
518 if(swap) SwapBytes((
char*)&one,
sizeof(
int), 1);
519 LOG(INFO) <<
"one after swap : " << one <<
"\n";
520 LOG(INFO) <<
"Swapping bytes from binary file (to be done)\n";
527 CHECK( std::string( __buf ) ==
"$EndMeshFormat" )
528 <<
"invalid file format entry "
530 <<
" instead of $EndMeshFormat\n";
532 DVLOG(2) <<
"[importergmsh] " << __buf <<
" (expect $PhysicalNames)\n";
534 std::vector<MeshMarkerName> meshMarkerNameMap = markerMap(MeshType::nDim);
535 if ( std::string( __buf ) ==
"$PhysicalNames" )
540 for (
int n = 0; n < nnames; ++n )
545 if ( boost::lexical_cast<double>( this->version() ) >= 2.1 )
547 __is >> topodim >>
id >> name;
548 DVLOG(2) <<
"[importergmsh] reading topodim: " << topodim <<
" id: " <<
id <<
" name: " << name <<
"\n";
551 else if ( this->version() ==
"2.0" )
555 boost::trim_if( name,boost::is_any_of(
"\"" ) );
557 if ( meshMarkerNameMap.empty() )
559 std::vector<int> data = {id, topodim};
560 mesh->addMarkerName( name,
id, topodim );
562 if ( M_ignorePhysicalName.find( name )!=M_ignorePhysicalName.end() ) this->setIgnorePhysicalGroup(
id );
564 if ( meshMarkerNameMap.empty() )
566 FEELPP_ASSERT( mesh->markerNames().size() == (
size_type )nnames )( mesh->markerNames().size() )( nnames ).error(
"invalid number of physical names" );
569 FEELPP_ASSERT( std::string( __buf ) ==
"$EndPhysicalNames" )
571 (
"$EndPhysicalNames" ).error (
"invalid file format entry" );
575 for(
auto it = meshMarkerNameMap.begin(), en = meshMarkerNameMap.end(); it != en; ++ it )
577 mesh->addMarkerName( it->name, it->ids[0], it->ids[1] );
584 DVLOG(2) <<
"buf: "<< __buf <<
"\n";
585 CHECK( std::string( __buf ) ==
"$NOD" ||
586 std::string( __buf ) ==
"$Nodes" ||
587 std::string( __buf ) ==
"$ParametricNodes" )
588 <<
"invalid nodes string '" << __buf <<
"' in gmsh importer. It should be either $NOD or $Nodes or $ParametricNodes.\n";
589 bool has_parametric_nodes = ( std::string( __buf ) ==
"$ParametricNodes" );
598 std::map<int, Feel::detail::GMSHPoint > gmshpts;
599 LOG(INFO) <<
"Reading "<< __n <<
" nodes\n";
603 for ( uint __i = 0; __i < __n; ++__i )
615 __is.read( (
char*)&
id,
sizeof(
int) );
616 if(swap) SwapBytes((
char*)&
id,
sizeof(
int), 1);
617 __is.read( (
char*)&x[0], 3*
sizeof(
double) );
618 if(swap) SwapBytes((
char*)&x[0],
sizeof(
double), 3);
624 if ( has_parametric_nodes )
627 gmshpts[id].parametric =
true;
628 CHECK( !binary ) <<
"GMSH Binary format not yet supported for parametric nodes\n";
629 int gdim = 0, gtag = 0;
632 __is >> gdim >> gtag;
640 else if ( gdim == 2 )
641 __is >> uv[0] >> uv[1];
643 gmshpts[id].gdim = gdim;
644 gmshpts[id].gtag = gtag;
659 DVLOG(2) <<
"buf: "<< __buf <<
"\n";
661 CHECK( std::string( __buf ) ==
"$ENDNOD" ||
662 std::string( __buf ) ==
"$EndNodes" ||
663 std::string( __buf ) ==
"$EndParametricNodes" )
664 <<
"invalid end nodes string '" << __buf
665 <<
"' in gmsh importer. It should be either $ENDNOD or $EndNodes or $EndParametricNodes\n";
671 CHECK( std::string( __buf ) ==
"$ELM" ||
672 std::string( __buf ) ==
"$Elements" )
673 <<
"invalid elements string " << __buf <<
" in gmsh importer, it should be either $ELM or $Elements\n";
682 LOG(INFO) <<
"Reading " << numElements <<
" elements...\n";
683 std::list<Feel::detail::GMSHElement> __et;
684 std::map<int,int> __idGmshToFeel;
685 std::map<int,int> __gt;
689 for(
int i = 0; i < numElements; i++)
691 int num, type, physical = 0, elementary = 0, parent = 0;
692 int dom1 = 0, dom2 = 0, numVertices;
693 std::vector<int> ghosts;
697 int partition = (this->worldComm().globalSize()>1)?this->worldComm().localRank():0;
702 int numPartitions = 1;
704 for(
int j = 0; j < numTags; j++)
708 if(j == 0) physical = tag;
709 else if(j == 1) elementary = tag;
710 else if(version < 2.2 && j == 2) partition = tag;
711 else if(version >= 2.2 && j == 2 && numTags > 3) numPartitions = tag;
712 else if(version >= 2.2 && j == 3) partition = tag-1;
713 else if(j >= 4 && j < 4 + numPartitions - 1) ghosts.push_back((-tag)-1);
714 else if(j == 3 + numPartitions && (numTags == 4 + numPartitions))
716 else if(j == 3 + numPartitions && (numTags == 5 + numPartitions)) {
721 CHECK(type != MSH_POLYG_ && type != MSH_POLYH_ && type != MSH_POLYG_B)
722 <<
"GMSH Element type " << type <<
" not supported by Feel++\n";
723 numVertices = MElement::getInfoMSH(type);
724 CHECK(numVertices!=0) <<
"Unknown number of vertices for element type " << type <<
"\n";
726 std::vector<int> indices(numVertices);
727 for(
int j = 0; j < numVertices; j++)
732 indices[j] = itoii[ indices[j] ];
735 if ( M_use_elementary_region_as_physical_region )
737 physical = elementary;
740 Feel::detail::GMSHElement gmshElt( num, type, physical, elementary,
741 numPartitions, partition, ghosts,
743 numVertices, indices,
744 this->worldComm().localRank(),this->worldComm().localSize() );
746 if ( gmshElt.isOnProcessor() ==
false ||
747 gmshElt.isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
750 __et.push_back( gmshElt );
752 if ( __gt.find( type ) != __gt.end() )
762 int numElementsPartial = 0;
763 while(numElementsPartial < numElements)
767 __is.read( (
char*)&header, 3*
sizeof(
int) );
768 if(swap) SwapBytes((
char*)header,
sizeof(
int), 3);
770 int type = header[0];
771 int numElems = header[1];
772 int numTags = header[2];
774 CHECK( type < MSH_NUM_TYPE ) <<
"Invalid GMSH element type " << type <<
"\n";
775 int numVertices = MElement::getInfoMSH(type,&name);
776 CHECK( numVertices > 0 ) <<
"Unsupported element type " << name <<
"\n";
778 unsigned int n = 1 + numTags + numVertices;
779 std::vector<int> data(n);
780 std::vector<int> indices( numVertices );
781 std::vector<int> ghosts;
783 for(
int i = 0; i < numElems; i++)
786 __is.read( (
char*)data.data(),
sizeof(int)*n );
787 if(swap) SwapBytes((
char*)data.data(),
sizeof(int), n);
791 int physical = (numTags > 0) ? data[1] : 0;
792 int elementary = (numTags > 1) ? data[2] : 0;
793 int numPartitions = (version >= 2.2 && numTags > 3) ? data[3] : 1;
794 int partition = (version < 2.2 && numTags > 2) ? data[3]-1 :
795 (version >= 2.2 && numTags > 3) ? data[4]-1 : 0;
796 if(numPartitions > 1)
797 for(
int j = 0; j < numPartitions - 1; j++)
798 ghosts.push_back( (-data[5 + j]) -1 );
799 int parent = (version < 2.2 && numTags > 3) ||
800 (version >= 2.2 && numPartitions && numTags > 3 + numPartitions) ||
801 (version >= 2.2 && !numPartitions && numTags > 2) ?
803 int dom1 = 0, dom2 = 0;
805 std::copy( &data[numTags + 1], &data[numTags + 1]+numVertices, indices.begin() );
809 for(
int j = 0; j < numVertices; j++)
811 indices[j] = itoii[ indices[j] ];
814 if ( M_use_elementary_region_as_physical_region )
816 physical = elementary;
819 Feel::detail::GMSHElement gmshElt( num, type, physical, elementary,
820 numPartitions, partition, ghosts,
822 numVertices, indices,
823 this->worldComm().localRank(),this->worldComm().localSize() );
825 if ( gmshElt.isOnProcessor() ==
false ||
826 gmshElt.isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
829 __et.push_back( gmshElt );
831 if ( __gt.find( type ) != __gt.end() )
837 numElementsPartial += numElems;
840 CHECK( numElementsPartial == numElements ) <<
"Invalid number of elements read from GMSH, read " << numElementsPartial <<
" element but expected " << numElements <<
"\n";
843 auto it = __gt.begin();
844 auto en = __gt.end();
845 for( ; it != en; ++it )
848 MElement::getInfoMSH( it->first, &name );
849 LOG(INFO) <<
"Read " << it->second <<
" " << name <<
" elements\n";
857 CHECK( std::string( __buf ) ==
"$ENDELM" ||
858 std::string( __buf ) ==
"$EndElements" )
859 <<
"invalid end elements string " << __buf
860 <<
" in gmsh importer. It should be either $ENDELM or $EndElements\n";
866 for ( uint __i = 0; __i < numElements; ++__i )
870 if ( __et[__i].isOnProcessor() ==
false ||
871 __et[__i].isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
874 switch ( __et[__i].type )
877 if ( mesh_type::nDim == 1 )
879 gmshpts[ __et[__i].indices[0] ].onbdy =
true;
889 if ( mesh_type::nDim == 2 )
891 for ( uint16_type jj = 0; jj < npoints_per_edge; ++jj )
893 gmshpts[ __et[__i].indices[jj] ].onbdy =
true;
906 if ( mesh_type::nDim == 3 )
908 for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
910 gmshpts[ __et[__i].indices[jj] ].onbdy =
true;
921 std::map<int,boost::tuple<int,int> > mapGhostElt;
923 node_type coords( mesh_type::nRealDim );
927 auto it_gmshElt = __et.begin();
928 auto const en_gmshElt = __et.end();
929 for ( ; it_gmshElt!=en_gmshElt ; ++it_gmshElt )
934 if ( it_gmshElt->isOnProcessor() ==
false ||
935 it_gmshElt->isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
938 for ( uint16_type p = 0; p < it_gmshElt->numVertices; ++p )
940 int ptid = it_gmshElt->indices[p];
942 if ( mesh->hasPoint( ptid ) )
945 auto const& gmshpt = gmshpts.find(ptid)->second;
946 for ( uint16_type j = 0; j < mesh_type::nRealDim; ++j )
947 coords[j] = gmshpt.x[j];
949 point_type pt( ptid, coords, gmshpt.onbdy );
951 if ( gmshpt.parametric )
953 pt.setGDim( gmshpt.gdim );
954 pt.setGTag( gmshpt.gtag );
956 if ( gmshpt.gdim < 3 )
958 pt.setParametricCoordinates( gmshpt.uv[0], gmshpt.uv[1] );
959 mesh->setParametric(
true );
962 mesh->addPoint( pt );
965 switch ( it_gmshElt->type )
970 addPoint( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
983 addEdge( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
997 addFace( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
1011 addVolume( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
1020 if ( it_gmshElt->isGhost() )
1022 mapGhostElt.insert( std::make_pair( it_gmshElt->num,boost::make_tuple( __idGmshToFeel[it_gmshElt->num], it_gmshElt->ghostPartitionId() ) ) );
1033 CHECK( mesh->numElements() > 0 ) <<
"The mesh does not have any elements.\n"
1034 <<
"something was not right with GMSH mesh importation.\n"
1035 <<
"please check that there are elements of topological dimension "
1036 << mesh_type::nDim <<
" in the mesh\n";
1038 if ( this->worldComm().localSize()>1 )
1039 updateGhostCellInfo( mesh, __idGmshToFeel, mapGhostElt );
1041 mesh->setNumVertices( std::accumulate( M_n_vertices.begin(), M_n_vertices.end(), 0,
1042 [](
int lhs, std::pair<int,int>
const& rhs )
1044 return lhs+rhs.second;
1046 if ( !mesh->markerNames().empty() &&
1047 ( mesh->markerNames().find(
"CrossPoints") != mesh->markerNames().end() ) &&
1048 ( mesh->markerNames().find(
"WireBasket") != mesh->markerNames().end() ) )
1054 DVLOG(2) <<
"done with reading and creating mesh from gmsh file\n";
1055 M_n_vertices.clear();
1058 template<
typename MeshType>
1062 addPoint( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1064 template<
typename MeshType>
1066 ImporterGmsh<MeshType>::addPoint( mesh_type*mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<1> )
1069 pf.setProcessIdInPartition( this->worldComm().localRank() );
1070 pf.setId( mesh->numFaces() );
1071 pf.setMarker( __e.physical );
1072 pf.setMarker2( __e.elementary );
1073 pf.setNumberOfPartitions( __e.numPartitions );
1074 pf.setProcessId( __e.partition );
1075 pf.setNeighborPartitionIds( __e.ghosts );
1077 pf.setPoint( 0, mesh->point( __e.indices[0] ) );
1079 if ( mesh->point( __e.indices[0] ).isOnBoundary() )
1081 pf.setOnBoundary(
true );
1083 M_n_vertices[ __e.indices[0] ] = 1;
1089 boost::tie( fit, inserted ) = mesh->addFace( pf );
1090 __idGmshToFeel=pf.id();
1092 DVLOG(2) <<
"added point on boundary ("
1093 << fit->isOnBoundary() <<
") with id :" << fit->id() <<
" and marker " << pf.marker()
1094 <<
" n1: " << mesh->point( __e.indices[0] ).node() <<
"\n";
1096 template<
typename MeshType>
1098 ImporterGmsh<MeshType>::addPoint( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<2> )
1100 auto pit = mesh->points().modify( mesh->pointIterator(__e.indices[0]),
1101 [&__e]( point_type& e )
1103 e.setMarker( __e.physical );
1104 e.setMarker2( __e.elementary );
1105 e.setNumberOfPartitions( __e.numPartitions );
1106 e.setProcessId( __e.partition );
1107 e.setNeighborPartitionIds( __e.ghosts );
1109 DVLOG(2) <<
"added point with id :" << mesh->pointIterator(__e.indices[0])->id() <<
" and marker " << mesh->pointIterator(__e.indices[0])->marker()
1110 <<
" n1: " << mesh->point( __e.indices[0] ).node() <<
"\n";
1113 template<
typename MeshType>
1115 ImporterGmsh<MeshType>::addPoint( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<3> )
1117 auto pit = mesh->pointIterator(__e.indices[0]);
1118 bool mod = mesh->points().modify( pit,
1119 [&__e]( point_type& e )
1121 e.setMarker( __e.physical );
1122 e.setMarker2( __e.elementary );
1123 e.setNumberOfPartitions( __e.numPartitions );
1124 e.setProcessId( __e.partition );
1125 e.setNeighborPartitionIds( __e.ghosts );
1127 DVLOG(2) <<
"added point (modified: " << mod <<
")with id :" << pit->id() <<
" and marker " << pit->marker()
1128 <<
" n1: " << mesh->point( __e.indices[0] ).node() <<
"\n";
1131 template<
typename MeshType>
1133 ImporterGmsh<MeshType>::addEdge( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel )
1135 addEdge( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1137 template<
typename MeshType>
1139 ImporterGmsh<MeshType>::addEdge( mesh_type*mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<1> )
1143 e.setProcessIdInPartition( this->worldComm().localRank() );
1144 e.setMarker( __e.physical );
1145 e.setMarker2( __e.elementary );
1146 e.setNumberOfPartitions( __e.numPartitions );
1147 e.setProcessId( __e.partition );
1148 e.setNeighborPartitionIds( __e.ghosts );
1156 int count_pt_on_boundary = 0;
1157 for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1159 e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1161 if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1162 ++count_pt_on_boundary;
1164 if ( count_pt_on_boundary >= 1 )
1166 e.setOnBoundary(
true );
1171 mesh->addElement( e );
1172 __idGmshToFeel=e.id();
1174 M_n_vertices[ __e.indices[0] ] = 1;
1175 M_n_vertices[ __e.indices[1] ] = 1;
1176 DVLOG(2) <<
"added edge with id :" << e.id()
1177 <<
" n1: " << mesh->point( __e.indices[0] ).node()
1178 <<
" n2: " << mesh->point( __e.indices[1] ).node() <<
"\n";
1181 template<
typename MeshType>
1183 ImporterGmsh<MeshType>::addEdge( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<2> )
1186 e.setProcessIdInPartition( this->worldComm().localRank() );
1187 e.setId( mesh->numFaces() );
1188 e.setMarker( __e.physical );
1189 e.setMarker2( __e.elementary );
1190 e.setNumberOfPartitions( __e.numPartitions );
1191 e.setProcessId( __e.partition );
1192 e.setNeighborPartitionIds( __e.ghosts );
1200 DCHECK( __e.indices.size() == npoints_per_edge )
1201 <<
"Invalid element indices, "
1202 <<
"indices size : " << __e.indices.size() <<
" points per edge : " << npoints_per_edge;
1203 int count_pt_on_boundary = 0;
1204 for ( uint16_type jj = 0; jj < npoints_per_edge; ++jj )
1206 e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1207 if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1208 ++count_pt_on_boundary;
1210 if ( count_pt_on_boundary >= 2 )
1211 e.setOnBoundary(
true );
1215 M_n_vertices[ __e.indices[0] ] = 1;
1216 M_n_vertices[ __e.indices[1] ] = 1;
1223 boost::tie( fit, inserted ) = mesh->addFace( e );
1224 __idGmshToFeel=e.id();
1226 DVLOG(2) <<
"added edge on boundary ("
1227 << fit->isOnBoundary() <<
") with id :" << fit->id()
1228 <<
" n1: " << mesh->point( __e.indices[0] ).node()
1229 <<
" n2: " << mesh->point( __e.indices[1] ).node() <<
"\n";
1231 template<
typename MeshType>
1233 ImporterGmsh<MeshType>::addEdge( mesh_type*mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<3> )
1236 e.setProcessIdInPartition( this->worldComm().localRank() );
1237 e.setId( mesh->numEdges() );
1238 e.setMarker( __e.physical );
1239 e.setMarker2( __e.elementary );
1240 e.setNumberOfPartitions( __e.numPartitions );
1241 e.setProcessId( __e.partition );
1242 e.setNeighborPartitionIds( __e.ghosts );
1250 int count_pt_on_boundary = 0;
1251 for ( uint16_type jj = 0; jj < npoints_per_edge; ++jj )
1253 e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1254 if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1255 ++count_pt_on_boundary;
1257 if ( count_pt_on_boundary >= 2 )
1259 e.setOnBoundary(
true );
1263 M_n_vertices[ __e.indices[0] ] = 1;
1264 M_n_vertices[ __e.indices[1] ] = 1;
1269 auto eit = mesh->addEdge( e );
1270 __idGmshToFeel=eit.id();
1272 if ( npoints_per_edge == 2 )
1273 DVLOG(2) <<
"added edge on boundary ("
1274 << eit.isOnBoundary() <<
") with id :" << eit.id()
1275 <<
" n1: " << eit.point( 0 ).node()
1276 <<
" n2: " << eit.point( 1 ).node() <<
"\n";
1280 template<
typename MeshType>
1282 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel )
1284 addFace( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1287 template<
typename MeshType>
1289 ImporterGmsh<MeshType>::addFace( mesh_type*, Feel::detail::GMSHElement
const&,
int & , mpl::int_<1> )
1292 template<
typename MeshType>
1294 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<2> )
1296 GmshOrdering<element_type> ordering;
1299 e.setProcessIdInPartition( this->worldComm().localRank() );
1300 e.setMarker( __e.physical );
1301 e.setMarker2( __e.elementary );
1302 e.setNumberOfPartitions( __e.numPartitions );
1303 e.setProcessId( __e.partition );
1304 e.setNeighborPartitionIds( __e.ghosts );
1314 int count_pt_on_boundary = 0;
1315 for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1318 if (!e.isGhostCell()) mesh->points().modify( mesh->pointIterator( __e.indices[jj] ), Feel::detail::UpdateProcessId(e.processId()) );
1319 e.setPoint( ordering.fromGmshId( jj ), mesh->point( __e.indices[jj] ) );
1320 if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1321 ++count_pt_on_boundary;
1328 count_pt_on_boundary >= 2 )
1330 e.setOnBoundary(
true );
1334 count_pt_on_boundary >= 2 )
1335 e.setOnBoundary(
true );
1339 mesh->addElement( e );
1340 __idGmshToFeel=e.id();
1342 M_n_vertices[ __e.indices[0] ] = 1;
1343 M_n_vertices[ __e.indices[1] ] = 1;
1344 M_n_vertices[ __e.indices[2] ] = 1;
1348 M_n_vertices[ __e.indices[3] ] = 1;
1350 template<
typename MeshType>
1352 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<3> )
1354 GmshOrdering<face_type> ordering;
1357 e.setProcessIdInPartition( this->worldComm().localRank() );
1358 e.setId( mesh->numFaces() );
1359 e.setMarker( __e.physical );
1360 e.setMarker2( __e.elementary );
1361 e.setNumberOfPartitions( __e.numPartitions );
1362 e.setProcessId( __e.partition );
1363 e.setNeighborPartitionIds( __e.ghosts );
1368 e.setOnBoundary(
true );
1378 for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
1381 e.setPoint( ordering.fromGmshId( jj ), mesh->point( __e.indices[jj] ) );
1389 boost::tie( fit, inserted ) = mesh->addFace( e );
1391 __idGmshToFeel=e.id();
1393 M_n_vertices[ __e.indices[0] ] = 1;
1394 M_n_vertices[ __e.indices[1] ] = 1;
1395 M_n_vertices[ __e.indices[2] ] = 1;
1399 M_n_vertices[ __e.indices[3] ] = 1;
1403 template<
typename MeshType>
1405 ImporterGmsh<MeshType>::addVolume( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel )
1407 addVolume( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1409 template<
typename MeshType>
1411 ImporterGmsh<MeshType>::addVolume( mesh_type*, Feel::detail::GMSHElement
const&,
int & , mpl::int_<1> )
1413 template<
typename MeshType>
1415 ImporterGmsh<MeshType>::addVolume( mesh_type*, Feel::detail::GMSHElement
const&,
int & , mpl::int_<2> )
1417 template<
typename MeshType>
1419 ImporterGmsh<MeshType>::addVolume( mesh_type* mesh, Feel::detail::GMSHElement
const& __e,
int & __idGmshToFeel, mpl::int_<3> )
1422 e.setProcessIdInPartition( this->worldComm().localRank() );
1423 GmshOrdering<element_type> ordering;
1424 e.setMarker( __e.physical );
1425 e.setMarker2( __e.elementary );
1426 e.setNumberOfPartitions( __e.numPartitions );
1427 e.setProcessId( __e.partition );
1428 e.setNeighborPartitionIds( __e.ghosts );
1442 int count_pt_on_boundary = 0;
1443 for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1446 if (!e.isGhostCell()) mesh->points().modify( mesh->pointIterator( __e.indices[jj] ), Feel::detail::UpdateProcessId(e.processId()) );
1448 e.setPoint( ordering.fromGmshId( jj ), mesh->point( __e.indices[jj] ) );
1450 if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1451 ++count_pt_on_boundary;
1459 && count_pt_on_boundary >= 3)
1460 e.setOnBoundary(
true );
1464 && count_pt_on_boundary >= 4)
1465 e.setOnBoundary(
true );
1468 mesh->addElement( e );
1469 __idGmshToFeel=e.id();
1471 M_n_vertices[ __e.indices[0] ] = 1;
1472 M_n_vertices[ __e.indices[1] ] = 1;
1473 M_n_vertices[ __e.indices[2] ] = 1;
1474 M_n_vertices[ __e.indices[3] ] = 1;
1478 M_n_vertices[ __e.indices[4] ] = 1;
1479 M_n_vertices[ __e.indices[5] ] = 1;
1480 M_n_vertices[ __e.indices[6] ] = 1;
1481 M_n_vertices[ __e.indices[7] ] = 1;
1485 template<
typename MeshType>
1487 ImporterGmsh<MeshType>::updateGhostCellInfo( mesh_type* mesh, std::map<int,int>
const& __idGmshToFeel, std::map<
int,boost::tuple<int,int> >
const& __mapGhostElt )
1490 std::vector<int> nbMsgToSend( this->worldComm().localSize() );
1491 std::fill( nbMsgToSend.begin(),nbMsgToSend.end(),0 );
1494 std::vector< std::map<int,int> > mapMsg( this->worldComm().localSize() );
1497 auto it_map = __mapGhostElt.begin();
1498 auto en_map = __mapGhostElt.end();
1500 for (
int cpt=0; it_map!=en_map; ++it_map,++cpt )
1502 auto idGmsh = it_map->first;
1503 auto idProc = it_map->second.template get<1>();
1505 std::cout <<
"[updateGhostCellInfo]----1---\n"
1506 <<
"I am the proc " << this->worldComm().globalRank()
1507 <<
" local proc " << this->worldComm().localRank()
1508 <<
" I send to the proc " << idProc <<
" for idGmsh " << idGmsh+1
1509 <<
" with tag "<< nbMsgToSend[idProc]
1510 <<
" the G " << mesh->element( it_map->second.template get<0>(),idProc ).G()
1514 this->worldComm().localComm().send( idProc , nbMsgToSend[idProc], idGmsh );
1516 mapMsg[idProc].insert( std::make_pair( nbMsgToSend[idProc],it_map->second.template get<0>() ) );
1518 nbMsgToSend[idProc]++;
1522 std::vector<int> nbMsgToRecv;
1523 mpi::all_to_all( this->worldComm().localComm(),
1528 for (
int proc=0; proc<this->worldComm().localSize(); ++proc )
1530 for (
int cpt=0 ; cpt<nbMsgToRecv[proc] ; ++cpt )
1534 this->worldComm().localComm().recv( proc, cpt, idGmsh );
1536 std::cout <<
"[updateGhostCellInfo]----2---\n"
1537 <<
"I am the proc" << this->worldComm().localRank()
1538 <<
" I receive of the proc " << proc <<
" for idGmsh " << idGmsh+1
1539 <<
" with tag "<< cpt
1540 <<
" the G " << mesh->element( __idGmshToFeel.find(idGmsh)->second ).G()
1541 <<
" with idFeel Classic " << __idGmshToFeel.find(idGmsh)->second
1545 this->worldComm().localComm().send( proc, cpt, __idGmshToFeel.find( idGmsh )->second );
1550 for (
int proc=0; proc<this->worldComm().localSize(); ++proc )
1552 for (
int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1556 this->worldComm().localComm().recv( proc, cpt, idFeel );
1558 auto elttt = mesh->elementIterator( mapMsg[proc][cpt],proc );
1559 mesh->elements().modify( elttt, Feel::detail::updateIdInOthersPartitions( proc, idFeel ) );
1561 std::cout <<
"[updateGhostCellInfo]----3---\n"
1562 <<
"END! I am the proc" << this->worldComm().localRank()<<
" I receive of the proc " << proc
1563 <<
" with tag "<< cpt
1564 <<
" for the G " << mesh->element( mapMsg[proc][cpt], proc ).G()
1565 <<
" with idFeel Classic " << mapMsg[proc][cpt]
1566 <<
" with idFeel " << idFeel
1567 <<
" and modif " << mesh->element( mapMsg[proc][cpt] , proc ).idInPartition( proc )