Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
importergmsh.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-11-16
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007,2008 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
31 #ifndef __ImporterGmsh_H
32 #define __ImporterGmsh_H 1
33 
34 #include <iostream>
35 #include <fstream>
36 #include <iomanip>
37 #include <algorithm>
38 
42 #include <boost/algorithm/string/trim.hpp>
43 
44 // Gmsh
45 #include <GModel.h>
46 #include <MElement.h>
47 
48 // from Gmsh
49 void SwapBytes(char *array, int size, int n);
50 
51 namespace Feel
52 {
53 namespace detail
54 {
55 class GMSHPoint
56 {
57 public:
58  Eigen::Vector3d x;
59  Eigen::Vector2d uv;
60  int id;
61  bool onbdy;
62  bool parametric;
63  int gdim,gtag;
64 
65  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
66 
67  GMSHPoint()
68  :
69  x( Eigen::Vector3d::Zero() ),
70  uv( Eigen::Vector2d::Zero() ),
71  id( 0 ),
72  onbdy( false ),
73  parametric( false ),
74  gdim( 0 ),
75  gtag( 0 )
76  {}
77  GMSHPoint(GMSHPoint const& p )
78  :
79  x( p.x ),
80  uv( p.uv ),
81  id ( p.id ),
82  onbdy( p.onbdy ),
83  parametric( p.parametric ),
84  gdim( p.gdim ),
85  gtag( p.gtag )
86  {}
87 };
88 struct GMSHElement
89 {
90  GMSHElement()
91  :
92  num( 0 ),
93  type( MSH_PNT ),
94  physical( 0 ),
95  elementary( 0 ),
96  numPartitions( 1 ),
97  partition( 0 ),
98  ghosts(),
99  is_on_processor( false ),
100  is_ghost( false ),
101  ghost_partition_id( -1 ),
102  parent( 0 ),
103  dom1(0), dom2(0),
104  numVertices(0),
105  indices()
106  {}
107 
108  GMSHElement( int n,
109  int t,
110  int p,
111  int e,
112  int _numPartitions,
113  int _partition,
114  std::vector<int> const& _ghosts,
115  int _parent,
116  int _dom1, int _dom2,
117  int _numVertices,
118  std::vector<int> const& _indices,
119  int worldcommrank,int worldcommsize)
120  :
121  num( n ),
122  type( t ),
123  physical( p ),
124  elementary( e ),
125  numPartitions( _numPartitions ),
126  partition( _partition ),
127  ghosts( _ghosts ),
128  is_on_processor( false ),
129  is_ghost( false ),
130  ghost_partition_id( partition ),
131  parent( _parent ),
132  dom1(_dom1), dom2(_dom2),
133  numVertices( _numVertices ),
134  indices( _indices )
135  {
136  if ( worldcommsize == 1 )
137  {
138  is_on_processor = true;
139  is_ghost = false;
140  }
141  else if ( worldcommrank == partition )
142  {
143  is_on_processor = true;
144  is_ghost = false;
145  }
146  else
147  {
148  // is the element a ghost cell
149  // look into ghosts if 'partition' is present
150  auto it = std::find( ghosts.begin(), ghosts.end(), worldcommrank );
151  if ( it != ghosts.end() )
152  {
153  is_on_processor = true;
154  is_ghost = true;
155  ghost_partition_id = partition;
156  }
157  }
158  }
159 
160  GMSHElement( GMSHElement const& g )
161  :
162  num( g.num ),
163  type( g.type ),
164  physical( g.physical ),
165  elementary( g.elementary ),
166  numPartitions( g.numPartitions ),
167  partition( g.partition ),
168  ghosts( g.ghosts ),
169  is_on_processor( g.is_on_processor ),
170  is_ghost( g.is_ghost ),
171  ghost_partition_id( g.ghost_partition_id ),
172  parent( g.parent ),
173  dom1( g.dom1 ), dom2( g.dom2 ),
174  numVertices( g.numVertices ),
175  indices( g.indices )
176  {}
177 
178  bool isOnProcessor() const { return is_on_processor; }
179  bool isGhost() const { return is_ghost; }
180  int ghostPartitionId() const { return ghost_partition_id; }
181 
182  template<typename IteratorBegin,typename IteratorEnd>
183  bool isIgnored(IteratorBegin it, IteratorEnd en ) const { return std::find( it, en, physical ) != en; }
184 
185  int num;
186  int type;
187  int physical;
188  int elementary;
189 
191  int numPartitions;
192  int partition;
193  std::vector<int> ghosts;
194  bool is_on_processor;
195  bool is_ghost;
196  int ghost_partition_id;
197 
198  int parent;
199  int dom1, dom2;
200 
201  // vertices
202  int numVertices;
203  std::vector<int> indices;
204 
205 
206 };
207 }
226 template<typename MeshType>
228  :
229 
230 public Importer<MeshType>
231 
232 {
233  typedef Importer<MeshType> super;
234 public:
235 
242  void setElementRegionAsPhysicalRegion( const bool param )
243  {
244  M_use_elementary_region_as_physical_region=param;
245  }
246 
250 
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;
257 
258  typedef typename mesh_type::face_iterator face_iterator;
260 
262 
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 );
267 
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 );
271 
272  static const uint16_type npoints_per_element = element_type::numPoints;
274 
278 
279  ImporterGmsh( WorldComm const& _worldcomm = Environment::worldComm() )
280  :
281  super( GMSH, _worldcomm ),
282  M_version( FEELPP_GMSH_FORMAT_VERSION ),
283  M_use_elementary_region_as_physical_region( false )
284  {
285  this->setIgnorePhysicalName( "FEELPP_GMSH_PHYSICALNAME_IGNORED" );
286  //showMe();
287  }
288 
289  explicit ImporterGmsh( std::string const& _fname, std::string _version = FEELPP_GMSH_FORMAT_VERSION,
290  WorldComm const& _worldcomm = Environment::worldComm() )
291  :
292  super( _fname, GMSH, _worldcomm ),
293  M_version( _version ),
294  M_use_elementary_region_as_physical_region( false )
295  {
296  this->setIgnorePhysicalName( "FEELPP_GMSH_PHYSICALNAME_IGNORED" );
297  //showMe();
298  }
299  ImporterGmsh( ImporterGmsh const & i )
300  :
301  super( 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 )
306  {
307  this->setIgnorePhysicalName( "FEELPP_GMSH_PHYSICALNAME_IGNORED" );
308  //showMe();
309  }
310  ~ImporterGmsh()
311  {}
312 
314 
318 
319 
321 
325 
329  std::string version() const
330  {
331  return M_version;
332  }
333 
338  boost::tuple<bool,boost::tuple<bool,int> > isElementOnProcessor( std::vector<int> const& tags ) const;
339 
341 
345 
346  void setVersion( std::string const& version )
347  {
348  M_version = version;
349  }
350 
351  void setIgnorePhysicalGroup( int i )
352  {
353  M_ignorePhysicalGroup.insert( i );
354  }
355  void setIgnorePhysicalName( std::string s )
356  {
357  M_ignorePhysicalName.insert( s );
358  }
359 
361 
365 
366  void visit( mesh_type* mesh );
367 
368  void showMe() const;
369 
371 
372 
373 
374 protected:
375 
376 private:
377 
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* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/, mpl::int_<2> );
381  void addPoint( mesh_type* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/, mpl::int_<3> );
382 
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* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/, mpl::int_<3> );
387 
388  void addFace( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel );
389  void addFace( mesh_type* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/, 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> );
392 
393  void addVolume( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel );
394  void addVolume( mesh_type* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/ , mpl::int_<1> );
395  void addVolume( mesh_type* /*mesh*/, Feel::detail::GMSHElement const& /*__e*/, int & /*__idGmshToFeel*/ , mpl::int_<2> );
396  void addVolume( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & /*__idGmshToFeel*/, mpl::int_<3> );
397 
398  void updateGhostCellInfo( mesh_type* mesh, std::map<int,int> const& __idGmshToFeel, std::map<int,boost::tuple<int,int> > const& __mapGhostElt );
399 
400 
401 private:
402 
403  std::string M_version;
404  std::map<int,int> M_n_vertices;
405  //std::vector<int> M_n_b_vertices;
406 
407  std::set<int> M_ignorePhysicalGroup;
408  std::set<std::string> M_ignorePhysicalName;
409  bool M_use_elementary_region_as_physical_region;
410 
411  //std::map<int,int> itoii;
412  //std::vector<int> ptseen;
413 
414 };
415 
416 
417 
418 template<typename MeshType>
419 void
420 ImporterGmsh<MeshType>::showMe() const
421 {
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";
425 }
426 template<typename MeshType>
427 boost::tuple<bool,boost::tuple<bool,int> >
428 ImporterGmsh<MeshType>::isElementOnProcessor( std::vector<int> const& tag ) const
429 {
430  bool is_element_on_processor = false;
431 
432  // if there is no partition tag (only 2 tags) then consider all elements on
433  // current processor
434  if ( tag.size() == 2 )
435  return boost::make_tuple( true,false );
436 
437  // tag[2] is the number of partition tags (1 in general 2 if cell share an
438  // face with 2 processors)
439  int rank = this->worldComm().localRank();
440  auto it = std::find_if( tag.begin()+3, tag.end(), [&rank] ( int i )
441  {
442  return i == rank;
443  } );
444 
445  if ( it != tag.end() )
446  is_element_on_processor = true;
447 
448  bool isGhostCell=false;
449  int idOfGhostCell=0;
450 
451  if ( tag[3]!=rank )
452  {
453  isGhostCell=true;
454  idOfGhostCell=tag[3];
455  }
456 
457  return boost::make_tuple( is_element_on_processor, boost::make_tuple( isGhostCell,idOfGhostCell ) );
458 }
459 template<typename MeshType>
460 void
462 {
463 
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" );
470 
471  DVLOG(2) << "visit(" << mesh_type::nDim << "D ) starts\n";
472  DVLOG(2) << "visit(" << mesh_type::nDim << "D ) filename = " << this->filename() << "\n";
473 
474  std::ifstream __is ( this->filename().c_str() );
475 
476  if ( !__is.is_open() )
477  {
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() );
482  }
483 
484  char __buf[256];
485  __is >> __buf;
486 
487  std::string theversion;
488  double version = 2.2;
489 
490  bool binary = false, swap = false;
491 
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" )
497  {
498 
499  // version file-type(0=ASCII,1=BINARY) data-size(sizeof(double))
500  int format, size;
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 );
506  if(format)
507  {
508  char c;
509  CHECK( (c=__is.get()) == '\n' ) << "Invalid character " << c << " should be newline\n";
510  binary = true;
511  LOG(INFO) << "GMSH mesh is in binary format\n";
512  int one;
513  __is.read( (char*)&one, sizeof(int) );
514  if(one != 1)
515  {
516  swap = true;
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";
521  }
522  //__is >> c;
523  //LOG(INFO) << "character= " << c << "\n";
524  }
525  // should be $EndMeshFormat
526  __is >> __buf;
527  CHECK( std::string( __buf ) == "$EndMeshFormat" )
528  << "invalid file format entry "
529  << __buf
530  << " instead of $EndMeshFormat\n";
531  __is >> __buf;
532  DVLOG(2) << "[importergmsh] " << __buf << " (expect $PhysicalNames)\n";
533 
534  std::vector<MeshMarkerName> meshMarkerNameMap = markerMap(MeshType::nDim);
535  if ( std::string( __buf ) == "$PhysicalNames" )
536  {
537  int nnames;
538  __is >> nnames;
539 
540  for ( int n = 0; n < nnames; ++n )
541  {
542  int id, topodim;
543  std::string name;
544 
545  if ( boost::lexical_cast<double>( this->version() ) >= 2.1 )
546  {
547  __is >> topodim >> id >> name;
548  DVLOG(2) << "[importergmsh] reading topodim: " << topodim << " id: " << id << " name: " << name << "\n";
549  }
550 
551  else if ( this->version() == "2.0" )
552  __is >> id >> name;
553 
554  boost::trim( name );
555  boost::trim_if( name,boost::is_any_of( "\"" ) );
556 
557  if ( meshMarkerNameMap.empty() )
558  {
559  std::vector<int> data = {id, topodim};
560  mesh->addMarkerName( name, id, topodim );
561  }
562  if ( M_ignorePhysicalName.find( name )!=M_ignorePhysicalName.end() ) this->setIgnorePhysicalGroup( id );
563  }
564  if ( meshMarkerNameMap.empty() )
565  {
566  FEELPP_ASSERT( mesh->markerNames().size() == ( size_type )nnames )( mesh->markerNames().size() )( nnames ).error( "invalid number of physical names" );
567  }
568  __is >> __buf;
569  FEELPP_ASSERT( std::string( __buf ) == "$EndPhysicalNames" )
570  ( __buf )
571  ( "$EndPhysicalNames" ).error ( "invalid file format entry" );
572  __is >> __buf;
573  }
574 
575  for(auto it = meshMarkerNameMap.begin(), en = meshMarkerNameMap.end(); it != en; ++ it )
576  {
577  mesh->addMarkerName( it->name, it->ids[0], it->ids[1] );
578  }
579  }
580 
581  //
582  // Read NODES
583  //
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" );
590  uint __n;
591  __is >> __n;
592 
593  // eat '\n' in binary mode otherwise the next binary read will get screwd
594  if ( binary )
595  __is.get();
596 
597 
598  std::map<int, Feel::detail::GMSHPoint > gmshpts;
599  LOG(INFO) << "Reading "<< __n << " nodes\n";
600 
601  Eigen::Vector3d x;
602  Eigen::Vector2d uv;
603  for ( uint __i = 0; __i < __n; ++__i )
604  {
605  int id = 0;
606  if ( !binary )
607  {
608  __is >> id
609  >> x[0]
610  >> x[1]
611  >> x[2];
612  }
613  else
614  {
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);
619 
620  }
621  gmshpts[id].id = id;
622  gmshpts[id].x = x;
623 
624  if ( has_parametric_nodes )
625  {
626 
627  gmshpts[id].parametric = true;
628  CHECK( !binary ) << "GMSH Binary format not yet supported for parametric nodes\n";
629  int gdim = 0, gtag = 0;
630  if ( !binary )
631  {
632  __is >> gdim >> gtag;
633 
634  // if gdim == 0 then u = v = 0
635  // if gdim == 3 then no need for a parametric point
636  // this logic is done later when filling the mesh data structure
637  if ( gdim == 1 )
638  __is >> uv[0];
639 
640  else if ( gdim == 2 )
641  __is >> uv[0] >> uv[1];
642  }
643  gmshpts[id].gdim = gdim;
644  gmshpts[id].gtag = gtag;
645  gmshpts[id].uv = uv;
646  }
647 
648  // stores mapping to be able to reorder the indices
649  // so that they are contiguous
650  //itoii[idpts[__i]] = __i;
651  }
652  //ptseen.resize( __n );
653  //std::fill( ptseen.begin(), ptseen.end(), -1 );
654  // eat '\n' in binary mode otherwise the next binary read will get screwd
655  if ( binary )
656  __is.get();
657 
658  __is >> __buf;
659  DVLOG(2) << "buf: "<< __buf << "\n";
660  // make sure that we have read all the points
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";
666 
667  //
668  // Read ELEMENTS
669  //
670  __is >> __buf;
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";
674 
675  int numElements;
676  __is >> numElements;
677 
678  // eat '\n' in binary mode otherwise the next binary read will get screwd
679  if ( binary )
680  __is.get();
681 
682  LOG(INFO) << "Reading " << numElements << " elements...\n";
683  std::list<Feel::detail::GMSHElement> __et; // tags in each element
684  std::map<int,int> __idGmshToFeel; // id Gmsh to id Feel
685  std::map<int,int> __gt;
686 
687  if ( !binary )
688  {
689  for(int i = 0; i < numElements; i++)
690  {
691  int num, type, physical = 0, elementary = 0, parent = 0;
692  int dom1 = 0, dom2 = 0, numVertices;
693  std::vector<int> ghosts;
694  int numTags;
695  // some faces may not be associated to a partition in the mesh file,
696  // hence will be read given the partition id 0 and will be discarded
697  int partition = (this->worldComm().globalSize()>1)?this->worldComm().localRank():0;
698  __is >> num // elm-number
699  >> type // elm-type
700  >> numTags; // number-of-tags
701 
702  int numPartitions = 1;
703 
704  for(int j = 0; j < numTags; j++)
705  {
706  int tag;
707  __is >> tag;
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))
715  parent = tag;
716  else if(j == 3 + numPartitions && (numTags == 5 + numPartitions)) {
717  dom1 = tag; j++;
718  __is >> dom2;
719  }
720  }
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";
725 
726  std::vector<int> indices(numVertices);
727  for(int j = 0; j < numVertices; j++)
728  {
729  __is >> indices[j];
730  // we do not renumber anymore
731 #if 0
732  indices[j] = itoii[ indices[j] ];
733 #endif
734  }
735  if ( M_use_elementary_region_as_physical_region )
736  {
737  physical = elementary;
738  }
739 
740  Feel::detail::GMSHElement gmshElt( num, type, physical, elementary,
741  numPartitions, partition, ghosts,
742  parent, dom1, dom2,
743  numVertices, indices,
744  this->worldComm().localRank(),this->worldComm().localSize() );
745 
746  if ( gmshElt.isOnProcessor() == false ||
747  gmshElt.isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
748  continue;
749 
750  __et.push_back( gmshElt );
751 
752  if ( __gt.find( type ) != __gt.end() )
753  ++__gt[ type ];
754  else
755  __gt[type]=1;
756 
757 
758  } // element description loop
759  } // !binary
760  else // binary case
761  {
762  int numElementsPartial = 0;
763  while(numElementsPartial < numElements)
764  {
765  int header[3];
766 
767  __is.read( (char*)&header, 3*sizeof(int) );
768  if(swap) SwapBytes((char*)header, sizeof(int), 3);
769 
770  int type = header[0];
771  int numElems = header[1];
772  int numTags = header[2];
773  char const* name;
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";
777 
778  unsigned int n = 1 + numTags + numVertices;
779  std::vector<int> data(n);
780  std::vector<int> indices( numVertices );
781  std::vector<int> ghosts;
782 
783  for(int i = 0; i < numElems; i++)
784  {
785  ghosts.clear();
786  __is.read( (char*)data.data(), sizeof(int)*n );
787  if(swap) SwapBytes((char*)data.data(), sizeof(int), n);
788 
789  //if(swap) SwapBytes((char*)data, sizeof(int), n);
790  int num = data[0];
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) ?
802  data[numTags] : 0;
803  int dom1 = 0, dom2 = 0;
804 
805  std::copy( &data[numTags + 1], &data[numTags + 1]+numVertices, indices.begin() );
806 
807  // we do not renumber anymore
808 #if 0
809  for(int j = 0; j < numVertices; j++)
810  {
811  indices[j] = itoii[ indices[j] ];
812  }
813 #endif
814  if ( M_use_elementary_region_as_physical_region )
815  {
816  physical = elementary;
817  }
818 
819  Feel::detail::GMSHElement gmshElt( num, type, physical, elementary,
820  numPartitions, partition, ghosts,
821  parent, dom1, dom2,
822  numVertices, indices,
823  this->worldComm().localRank(),this->worldComm().localSize() );
824 
825  if ( gmshElt.isOnProcessor() == false ||
826  gmshElt.isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
827  continue;
828 
829  __et.push_back( gmshElt );
830 
831  if ( __gt.find( type ) != __gt.end() )
832  ++__gt[ type ];
833  else
834  __gt[type]=1;
835 
836  }
837  numElementsPartial += numElems;
838 
839  } // while
840  CHECK( numElementsPartial == numElements ) << "Invalid number of elements read from GMSH, read " << numElementsPartial << " element but expected " << numElements << "\n";
841  } // binary
842 
843  auto it = __gt.begin();
844  auto en = __gt.end();
845  for( ; it != en; ++it )
846  {
847  const char* name;
848  MElement::getInfoMSH( it->first, &name );
849  LOG(INFO) << "Read " << it->second << " " << name << " elements\n";
850  }
851 
852  if ( binary )
853  __is >> __buf;
854 
855  // make sure that we have read everything
856  __is >> __buf;
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";
861 
862 #if 0
863  //
864  // FILL Mesh Data Structure
865  //
866  for ( uint __i = 0; __i < numElements; ++__i )
867  {
868  // if the element is not associated to the processor (in partition or ghost) or
869  // if the physical entity is ignored
870  if ( __et[__i].isOnProcessor() == false ||
871  __et[__i].isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
872  continue;
873 
874  switch ( __et[__i].type )
875  {
876  case GMSH_POINT:
877  if ( mesh_type::nDim == 1 )
878  {
879  gmshpts[ __et[__i].indices[0] ].onbdy = true;
880  }
881 
882  break;
883 
884  case GMSH_LINE:
885  case GMSH_LINE_2:
886  case GMSH_LINE_3:
887  case GMSH_LINE_4:
888  case GMSH_LINE_5:
889  if ( mesh_type::nDim == 2 )
890  {
891  for ( uint16_type jj = 0; jj < npoints_per_edge; ++jj )
892  {
893  gmshpts[ __et[__i].indices[jj] ].onbdy = true;
894  }
895  }
896 
897  break;
898 
899  case GMSH_QUADRANGLE:
900  case GMSH_QUADRANGLE_2:
901  case GMSH_TRIANGLE:
902  case GMSH_TRIANGLE_2:
903  case GMSH_TRIANGLE_3:
904  case GMSH_TRIANGLE_4:
905  case GMSH_TRIANGLE_5:
906  if ( mesh_type::nDim == 3 )
907  {
908  for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
909  {
910  gmshpts[ __et[__i].indices[jj] ].onbdy = true;
911  }
912  }
913 
914  break;
915 
916  default:
917  break;
918  }
919  }
920 #endif
921  std::map<int,boost::tuple<int,int> > mapGhostElt;
922 
923  node_type coords( mesh_type::nRealDim );
924 
925  //M_n_b_vertices.resize( __n );
926  //M_n_b_vertices.assign( __n, 0 );
927  auto it_gmshElt = __et.begin();
928  auto const en_gmshElt = __et.end();
929  for ( ; it_gmshElt!=en_gmshElt ; ++it_gmshElt )
930  // add the element to the mesh
931  {
932  // if the element is not associated to the processor (in partition or ghost) or
933  // if the physical entity is ignored
934  if ( it_gmshElt->isOnProcessor() == false ||
935  it_gmshElt->isIgnored(M_ignorePhysicalGroup.begin(), M_ignorePhysicalGroup.end()) )
936  continue;
937  // add the points associates to the element on the processor
938  for ( uint16_type p = 0; p < it_gmshElt->numVertices; ++p )
939  {
940  int ptid = it_gmshElt->indices[p];
941  // don't do anything if the point is already registered
942  if ( mesh->hasPoint( ptid ) )
943  continue;
944 
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];
948 
949  point_type pt( ptid, coords, gmshpt.onbdy );
950 
951  if ( gmshpt.parametric )
952  {
953  pt.setGDim( gmshpt.gdim );
954  pt.setGTag( gmshpt.gtag );
955 
956  if ( gmshpt.gdim < 3 )
957  {
958  pt.setParametricCoordinates( gmshpt.uv[0], gmshpt.uv[1] );
959  mesh->setParametric( true );
960  }
961  }
962  mesh->addPoint( pt );
963  } // loop over local points
964 
965  switch ( it_gmshElt->type )
966  {
967  // Points
968  case GMSH_POINT:
969  {
970  addPoint( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
971 
972 
973  break;
974  }
975 
976  // Edges
977  case GMSH_LINE:
978  case GMSH_LINE_2:
979  case GMSH_LINE_3:
980  case GMSH_LINE_4:
981  case GMSH_LINE_5:
982  {
983  addEdge( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
984 
985  break;
986  }
987 
988  // Faces
989  case GMSH_TRIANGLE:
990  case GMSH_TRIANGLE_2:
991  case GMSH_TRIANGLE_3:
992  case GMSH_TRIANGLE_4:
993  case GMSH_TRIANGLE_5:
994  case GMSH_QUADRANGLE:
995  case GMSH_QUADRANGLE_2:
996  {
997  addFace( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
998 
999  break;
1000  }
1001 
1002  // Volumes
1003  case GMSH_TETRAHEDRON:
1004  case GMSH_TETRAHEDRON_2:
1005  case GMSH_TETRAHEDRON_3:
1006  case GMSH_TETRAHEDRON_4:
1007  case GMSH_TETRAHEDRON_5:
1008  case GMSH_HEXAHEDRON:
1009  case GMSH_HEXAHEDRON_2:
1010  {
1011  addVolume( mesh, *it_gmshElt, __idGmshToFeel[it_gmshElt->num] );
1012 
1013  break;
1014  }
1015 
1016  default:
1017  break;
1018  }
1019 
1020  if ( it_gmshElt->isGhost() )
1021  {
1022  mapGhostElt.insert( std::make_pair( it_gmshElt->num,boost::make_tuple( __idGmshToFeel[it_gmshElt->num], it_gmshElt->ghostPartitionId() ) ) );
1023  }
1024 
1025  } // loop over geometric entities in gmsh file (can be elements or faces)
1026 
1027  if (VLOG_IS_ON(4))
1028  {
1029  //for(int i = 0; i < ptseen.size(); ++i )
1030  //if ( ptseen[i] == -1 )
1031  //LOG(WARNING) << "Point with id " << i << " not in element connectivity";
1032  }
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";
1037 
1038  if ( this->worldComm().localSize()>1 )
1039  updateGhostCellInfo( mesh, __idGmshToFeel, mapGhostElt );
1040 
1041  mesh->setNumVertices( std::accumulate( M_n_vertices.begin(), M_n_vertices.end(), 0,
1042  []( int lhs, std::pair<int,int> const& rhs )
1043  {
1044  return lhs+rhs.second;
1045  } ) );
1046  if ( !mesh->markerNames().empty() &&
1047  ( mesh->markerNames().find("CrossPoints") != mesh->markerNames().end() ) &&
1048  ( mesh->markerNames().find("WireBasket") != mesh->markerNames().end() ) )
1049  {
1050  //LOG(INFO) << "[substructuring] marker cp" << mesh->markerName("CrossPoints") << "\n";
1051  //LOG(INFO) << "[substructuring] marker wb" << mesh->markerName("WireBasket") << "\n";
1052  //LOG(INFO) << "[substructuring] n cp: " << std::distance( mesh->beginPointWithMarker( mesh->markerName("CrossPoints") ), mesh->endPointWithMarker( mesh->markerName("CrossPoints") ) ) << "\n";
1053  }
1054  DVLOG(2) << "done with reading and creating mesh from gmsh file\n";
1055  M_n_vertices.clear();
1056 }
1057 
1058 template<typename MeshType>
1059 void
1060 ImporterGmsh<MeshType>::addPoint( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel )
1061 {
1062  addPoint( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1063 }
1064 template<typename MeshType>
1065 void
1066 ImporterGmsh<MeshType>::addPoint( mesh_type*mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<1> )
1067 {
1068  face_type pf;
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 );
1076 
1077  pf.setPoint( 0, mesh->point( __e.indices[0] ) );
1078 // ptseen[mesh->point( __e.indices[0] ).id()]=1;
1079  if ( mesh->point( __e.indices[0] ).isOnBoundary() )
1080  {
1081  pf.setOnBoundary( true );
1082  }
1083  M_n_vertices[ __e.indices[0] ] = 1;
1084 
1085  //M_n_b_vertices[ __e.indices[0] ] = 1;
1086 
1087  face_iterator fit;
1088  bool inserted;
1089  boost::tie( fit, inserted ) = mesh->addFace( pf );
1090  __idGmshToFeel=pf.id();
1091 
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";
1095 }
1096 template<typename MeshType>
1097 void
1098 ImporterGmsh<MeshType>::addPoint( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<2> )
1099 {
1100  auto pit = mesh->points().modify( mesh->pointIterator(__e.indices[0]),
1101  [&__e]( point_type& e )
1102  {
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 );
1108  } );
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";
1111 
1112 }
1113 template<typename MeshType>
1114 void
1115 ImporterGmsh<MeshType>::addPoint( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<3> )
1116 {
1117  auto pit = mesh->pointIterator(__e.indices[0]);
1118  bool mod = mesh->points().modify( pit,
1119  [&__e]( point_type& e )
1120  {
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 );
1126  } );
1127  DVLOG(2) << "added point (modified: " << mod << ")with id :" << pit->id() << " and marker " << pit->marker()
1128  << " n1: " << mesh->point( __e.indices[0] ).node() << "\n";
1129 }
1130 
1131 template<typename MeshType>
1132 void
1133 ImporterGmsh<MeshType>::addEdge( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel )
1134 {
1135  addEdge( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1136 }
1137 template<typename MeshType>
1138 void
1139 ImporterGmsh<MeshType>::addEdge( mesh_type*mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<1> )
1140 {
1141  element_type e;
1142  //e.setWorldComm(this->worldComm());
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 );
1149 
1150  if ( __e.type == GMSH_LINE ||
1151  __e.type == GMSH_LINE_2 ||
1152  __e.type == GMSH_LINE_3 ||
1153  __e.type == GMSH_LINE_4 ||
1154  __e.type == GMSH_LINE_5 )
1155  {
1156  int count_pt_on_boundary = 0;
1157  for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1158  {
1159  e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1160  //ptseen[mesh->point( __e.indices[jj] ).id()]=1;
1161  if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1162  ++count_pt_on_boundary;
1163  }
1164  if ( count_pt_on_boundary >= 1 )
1165  {
1166  e.setOnBoundary( true );
1167  }
1168 
1169  }
1170 
1171  mesh->addElement( e );
1172  __idGmshToFeel=e.id();
1173 
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";
1179 }
1180 
1181 template<typename MeshType>
1182 void
1183 ImporterGmsh<MeshType>::addEdge( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<2> )
1184 {
1185  face_type e;
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 );
1193 
1194  if ( __e.type == GMSH_LINE ||
1195  __e.type == GMSH_LINE_2 ||
1196  __e.type == GMSH_LINE_3 ||
1197  __e.type == GMSH_LINE_4 ||
1198  __e.type == GMSH_LINE_5 )
1199  {
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 )
1205  {
1206  e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1207  if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1208  ++count_pt_on_boundary;
1209  }
1210  if ( count_pt_on_boundary >= 2 )
1211  e.setOnBoundary( true );
1212 
1213  }
1214 
1215  M_n_vertices[ __e.indices[0] ] = 1;
1216  M_n_vertices[ __e.indices[1] ] = 1;
1217 
1218  //M_n_b_vertices[ __e.indices[0] ] = 1;
1219  //M_n_b_vertices[ __e.indices[1] ] = 1;
1220 
1221  bool inserted;
1222  face_iterator fit;
1223  boost::tie( fit, inserted ) = mesh->addFace( e );
1224  __idGmshToFeel=e.id();
1225 
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";
1230 }
1231 template<typename MeshType>
1232 void
1233 ImporterGmsh<MeshType>::addEdge( mesh_type*mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<3> )
1234 {
1235  edge_type e;
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 );
1243 
1244  if ( __e.type == GMSH_LINE ||
1245  __e.type == GMSH_LINE_2 ||
1246  __e.type == GMSH_LINE_3 ||
1247  __e.type == GMSH_LINE_4 ||
1248  __e.type == GMSH_LINE_5 )
1249  {
1250  int count_pt_on_boundary = 0;
1251  for ( uint16_type jj = 0; jj < npoints_per_edge; ++jj )
1252  {
1253  e.setPoint( jj, mesh->point( __e.indices[jj] ) );
1254  if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1255  ++count_pt_on_boundary;
1256  }
1257  if ( count_pt_on_boundary >= 2 )
1258  {
1259  e.setOnBoundary( true );
1260  }
1261  }
1262 
1263  M_n_vertices[ __e.indices[0] ] = 1;
1264  M_n_vertices[ __e.indices[1] ] = 1;
1265 
1266  //M_n_b_vertices[ __e.indices[0] ] = 1;
1267  //M_n_b_vertices[ __e.indices[1] ] = 1;
1268 
1269  auto eit = mesh->addEdge( e );
1270  __idGmshToFeel=eit.id();
1271 
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";
1277 
1278 }
1279 
1280 template<typename MeshType>
1281 void
1282 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel )
1283 {
1284  addFace( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1285 }
1286 
1287 template<typename MeshType>
1288 void
1289 ImporterGmsh<MeshType>::addFace( mesh_type*, Feel::detail::GMSHElement const&, int & /*__idGmshToFeel*/, mpl::int_<1> )
1290 {}
1291 
1292 template<typename MeshType>
1293 void
1294 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<2> )
1295 {
1296  GmshOrdering<element_type> ordering;
1297 
1298  element_type e;
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 );
1305 
1306  if ( __e.type == GMSH_QUADRANGLE ||
1307  __e.type == GMSH_QUADRANGLE_2 ||
1308  __e.type == GMSH_TRIANGLE ||
1309  __e.type == GMSH_TRIANGLE_2 ||
1310  __e.type == GMSH_TRIANGLE_3 ||
1311  __e.type == GMSH_TRIANGLE_4 ||
1312  __e.type == GMSH_TRIANGLE_5 )
1313  {
1314  int count_pt_on_boundary = 0;
1315  for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1316  {
1317  //ptseen[mesh->point( __e.indices[jj] ).id()]=1;
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;
1322  }
1323  if ( ( __e.type == GMSH_TRIANGLE ||
1324  __e.type == GMSH_TRIANGLE_2 ||
1325  __e.type == GMSH_TRIANGLE_3 ||
1326  __e.type == GMSH_TRIANGLE_4 ||
1327  __e.type == GMSH_TRIANGLE_5 ) &&
1328  count_pt_on_boundary >= 2 )
1329  {
1330  e.setOnBoundary( true );
1331  }
1332  if ( (__e.type == GMSH_QUADRANGLE ||
1333  __e.type == GMSH_QUADRANGLE_2) &&
1334  count_pt_on_boundary >= 2 )
1335  e.setOnBoundary( true );
1336 
1337  }
1338 
1339  mesh->addElement( e );
1340  __idGmshToFeel=e.id();
1341 
1342  M_n_vertices[ __e.indices[0] ] = 1;
1343  M_n_vertices[ __e.indices[1] ] = 1;
1344  M_n_vertices[ __e.indices[2] ] = 1;
1345 
1346  if ( __e.type == GMSH_QUADRANGLE ||
1347  __e.type == GMSH_QUADRANGLE_2 )
1348  M_n_vertices[ __e.indices[3] ] = 1;
1349 }
1350 template<typename MeshType>
1351 void
1352 ImporterGmsh<MeshType>::addFace( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<3> )
1353 {
1354  GmshOrdering<face_type> ordering;
1355 
1356  face_type e;
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 );
1364 
1365  // we consider in 3D that all faces provided by Gmsh are on the boundary
1366  // which might not be true
1367  // \warning there might be a bug here
1368  e.setOnBoundary( true );
1369 
1370  if ( __e.type == GMSH_QUADRANGLE ||
1371  __e.type == GMSH_QUADRANGLE_2 ||
1372  __e.type == GMSH_TRIANGLE ||
1373  __e.type == GMSH_TRIANGLE_2 ||
1374  __e.type == GMSH_TRIANGLE_3 ||
1375  __e.type == GMSH_TRIANGLE_4 ||
1376  __e.type == GMSH_TRIANGLE_5 )
1377  {
1378  for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
1379  {
1380  //ptseen[mesh->point( __e.indices[jj] ).id()]=1;
1381  e.setPoint( ordering.fromGmshId( jj ), mesh->point( __e.indices[jj] ) );
1382  }
1383 
1384  //e.setPoint( jj, mesh->point( __e[jj] ) );
1385  }
1386 
1387  bool inserted;
1388  face_iterator fit;
1389  boost::tie( fit, inserted ) = mesh->addFace( e );
1390 
1391  __idGmshToFeel=e.id();
1392 
1393  M_n_vertices[ __e.indices[0] ] = 1;
1394  M_n_vertices[ __e.indices[1] ] = 1;
1395  M_n_vertices[ __e.indices[2] ] = 1;
1396 
1397  if ( __e.type == GMSH_QUADRANGLE ||
1398  __e.type == GMSH_QUADRANGLE_2 )
1399  M_n_vertices[ __e.indices[3] ] = 1;
1400 
1401 }
1402 
1403 template<typename MeshType>
1404 void
1405 ImporterGmsh<MeshType>::addVolume( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel )
1406 {
1407  addVolume( mesh, __e, __idGmshToFeel, mpl::int_<mesh_type::nDim>() );
1408 }
1409 template<typename MeshType>
1410 void
1411 ImporterGmsh<MeshType>::addVolume( mesh_type*, Feel::detail::GMSHElement const&, int & /*__idGmshToFeel*/, mpl::int_<1> )
1412 {}
1413 template<typename MeshType>
1414 void
1415 ImporterGmsh<MeshType>::addVolume( mesh_type*, Feel::detail::GMSHElement const&, int & /*__idGmshToFeel*/, mpl::int_<2> )
1416 {}
1417 template<typename MeshType>
1418 void
1419 ImporterGmsh<MeshType>::addVolume( mesh_type* mesh, Feel::detail::GMSHElement const& __e, int & __idGmshToFeel, mpl::int_<3> )
1420 {
1421  element_type e;
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 );
1429 
1430 
1431  //
1432  // WARNING: not yet done for high order in 3D !!!
1433  //
1434  if ( __e.type == GMSH_HEXAHEDRON ||
1435  __e.type == GMSH_HEXAHEDRON_2 ||
1436  __e.type == GMSH_TETRAHEDRON ||
1437  __e.type == GMSH_TETRAHEDRON_2 ||
1438  __e.type == GMSH_TETRAHEDRON_3 ||
1439  __e.type == GMSH_TETRAHEDRON_4 ||
1440  __e.type == GMSH_TETRAHEDRON_5 )
1441  {
1442  int count_pt_on_boundary = 0;
1443  for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1444  {
1445  //ptseen[mesh->point( __e.indices[jj] ).id()]=1;
1446  if (!e.isGhostCell()) mesh->points().modify( mesh->pointIterator( __e.indices[jj] ), Feel::detail::UpdateProcessId(e.processId()) );
1447  //std::cout << "gmsh index " << jj << " -> " << ordering.fromGmshId(jj) << " -> " << mesh->point( __e[jj] ).id()+1 << " : " << mesh->point( __e[jj] ).node() << "\n";
1448  e.setPoint( ordering.fromGmshId( jj ), mesh->point( __e.indices[jj] ) );
1449 
1450  if ( mesh->point( __e.indices[jj] ).isOnBoundary() )
1451  ++count_pt_on_boundary;
1452  }
1453  // the tet share a face with the boundary
1454  if ( ( __e.type == GMSH_TETRAHEDRON ||
1455  __e.type == GMSH_TETRAHEDRON_2 ||
1456  __e.type == GMSH_TETRAHEDRON_3 ||
1457  __e.type == GMSH_TETRAHEDRON_4 ||
1458  __e.type == GMSH_TETRAHEDRON_5 )
1459  && count_pt_on_boundary >= 3)
1460  e.setOnBoundary( true );
1461  // the hex share a face with the boundary
1462  if ( ( __e.type == GMSH_HEXAHEDRON ||
1463  __e.type == GMSH_HEXAHEDRON_2 )
1464  && count_pt_on_boundary >= 4)
1465  e.setOnBoundary( true );
1466  }
1467 
1468  mesh->addElement( e );
1469  __idGmshToFeel=e.id();
1470 
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;
1475 
1476  if ( __e.type == GMSH_HEXAHEDRON )
1477  {
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;
1482  }
1483 }
1484 
1485 template<typename MeshType>
1486 void
1487 ImporterGmsh<MeshType>::updateGhostCellInfo( mesh_type* mesh, std::map<int,int> const& __idGmshToFeel, std::map<int,boost::tuple<int,int> > const& __mapGhostElt )
1488 {
1489  // counter of msg sent for each process
1490  std::vector<int> nbMsgToSend( this->worldComm().localSize() );
1491  std::fill( nbMsgToSend.begin(),nbMsgToSend.end(),0 );
1492 
1493  // map usefull to get final result
1494  std::vector< std::map<int,int> > mapMsg( this->worldComm().localSize() );
1495 
1496  // iterate over ghost elt
1497  auto it_map = __mapGhostElt.begin();
1498  auto en_map = __mapGhostElt.end();
1499 
1500  for ( int cpt=0; it_map!=en_map; ++it_map,++cpt )
1501  {
1502  auto idGmsh = it_map->first;
1503  auto idProc = it_map->second.template get<1>();
1504 #if 0
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()
1511  << std::endl;
1512 #endif
1513  // send
1514  this->worldComm().localComm().send( idProc , nbMsgToSend[idProc], idGmsh );
1515  // save tag of request
1516  mapMsg[idProc].insert( std::make_pair( nbMsgToSend[idProc],it_map->second.template get<0>() ) );
1517  // update nb send
1518  nbMsgToSend[idProc]++;
1519  }
1520 
1521  // counter of msg received for each process
1522  std::vector<int> nbMsgToRecv;
1523  mpi::all_to_all( this->worldComm().localComm(),
1524  nbMsgToSend,
1525  nbMsgToRecv );
1526 
1527  // get gmsh id asked and re-send the correspond id Feel
1528  for ( int proc=0; proc<this->worldComm().localSize(); ++proc )
1529  {
1530  for ( int cpt=0 ; cpt<nbMsgToRecv[proc] ; ++cpt )
1531  {
1532  int idGmsh;
1533  //reception idGmsh
1534  this->worldComm().localComm().recv( proc, cpt, idGmsh );
1535 #if 0
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
1542  << std::endl;
1543 #endif
1544  //re-send idFeel
1545  this->worldComm().localComm().send( proc, cpt, __idGmshToFeel.find( idGmsh )->second );
1546  }
1547  }
1548 
1549  // get response to initial request and update Feel::Mesh data
1550  for ( int proc=0; proc<this->worldComm().localSize(); ++proc )
1551  {
1552  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1553  {
1554  int idFeel;
1555  // receive idFeel
1556  this->worldComm().localComm().recv( proc, cpt, idFeel );
1557  // update data
1558  auto elttt = mesh->elementIterator( mapMsg[proc][cpt],proc );
1559  mesh->elements().modify( elttt, Feel::detail::updateIdInOthersPartitions( proc, idFeel ) );
1560 #if 0
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 )
1568  << std::endl;
1569 #endif
1570  }
1571  }
1572 
1573 }
1574 
1575 } // Feel
1576 
1577 
1578 
1579 #endif /* __ImporterGmsh_H */

Generated on Sun Oct 20 2013 08:25:00 for Feel++ by doxygen 1.8.4