LifeV
ImporterMesh3D.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief Mesh reader from mesh3d files
30 
31  @author Luca Formaggia <luca.formaggia@polimi.it>
32  @contributor JFG, Nur Aiman Fadel <nur.fadel@mail.polimi.it>
33  @maintainer Nur Aiman Fadel <nur.fadel@mail.polimi.it>
34 
35  @date 29-06-2002
36 
37  Mesh reader that it is able to read 3d meshes.<br>
38  INRIAMesh used either spaces or CR as separators.<br>
39  */
40 
41 #ifndef _IMPORTERMESH3D_HH_
42 #define _IMPORTERMESH3D_HH_ 1
43 
44 
45 #include <boost/lambda/bind.hpp>
46 #include <boost/lambda/if.hpp>
47 #include <boost/lambda/lambda.hpp>
48 
49 
50 #include <lifev/core/LifeV.hpp>
51 
52 #include <lifev/core/util/StringUtility.hpp>
53 
54 #include <lifev/core/mesh/MeshElementBare.hpp>
55 
56 #include <lifev/core/mesh/MeshChecks.hpp>
57 #include <lifev/core/mesh/InternalEntitySelector.hpp>
58 #include <lifev/core/mesh/BareMesh.hpp>
59 #include <lifev/core/mesh/RegionMesh.hpp>
60 
61 namespace LifeV
62 {
63 
64 // @name Public typedefs
65 //@{
66 typedef boost::numeric::ublas::vector<Real> Vector;
67 //@}
68 
69 namespace
70 {
71 
72 struct FiveNumbers
73 {
74 public:
75  UInt i1, i2, i3, i4;
76  Int ibc;
77 };
78 
79 struct FaceHelp
80 {
81 public:
82  UInt i[9];
83  Int ibc;
84 };
85 
86 }
87 
88 // ===================================================
89 // Mpp mesh readers
90 // ===================================================
91 
92 //! readMppFileHead - reads mesh++ Tetra meshes.
93 /*!
94  It converts Tetra meshes into Quadratic Tetra if needed it.
95 
96  @param myStream,
97  @param numberVertices,
98  @param numberBoundaryVertices,
99  @param numberBoundaryFaces,
100  @param numberBoundaryEdges,
101  @param numberVolumes,
102  @return false if mesh check is unsuccessfull.
103 */
104 
105 bool
106 readMppFileHead ( std::ifstream& myStream,
107  UInt& numberVertices,
108  UInt& numberBoundaryVertices,
109  UInt& numberBoundaryFaces,
110  UInt& numberBoundaryEdges,
111  UInt& numberVolumes );
112 
113 //! readMppFile - reads mesh++ Tetra meshes.
114 /*!
115  It converts Tetra meshes into Quadratic Tetra if needed it.
116 
117  @param mesh, the mesh data structure to fill in.
118  @param fileName, the name of the mesh file to read.
119  @param regionFlag, the identifier for the region.
120  @param verbose, setting it as true, the output is verbose (the default is false).
121  @return true if everything went fine, false otherwise.
122 */
123 
124 template <typename GeoShape, typename MC>
125 bool
126 readMppFile ( RegionMesh<GeoShape, MC>& mesh,
127  const std::string& fileName,
128  markerID_Type regionFlag,
129  bool verbose = false )
130 {
131  typedef RegionMesh<GeoShape, MC> mesh_Type;
132 
133  const int idOffset = 1; //IDs in MPP files start from 1
134 
135  std::string line;
136 
137  Real x, y, z;
138 
139  Int ity, ity_id;
140 
141  UInt done = 0;
142  UInt i;
143  UInt numberVertices ( 0 ), numberBoundaryVertices ( 0 ),
144  numberFaces ( 0 ), numberBoundaryFaces ( 0 ),
145  numberPoints ( 0 ), numberBoundaryPoints ( 0 ),
146  numberEdges ( 0 ), numberBoundaryEdges ( 0 );
147  UInt numberVolumes ( 0 );
148  UInt p1, p2, p3, p4;
149 
150  std::stringstream discardedLog;
151 
152  std::ostream& oStr = verbose ? std::cout : discardedLog;
153 
154  ASSERT_PRE0 ( GeoShape::S_shape == TETRA, "Sorry, readMppFiles reads only tetra meshes" );
155  ASSERT_PRE0 ( GeoShape::S_numVertices <= 6, "Sorry, readMppFiles handles only liner&quad tetras" );
156 
157  // open stream to read header
158 
159  std::ifstream hstream ( fileName.c_str() );
160 
161  if ( hstream.fail() )
162  {
163  std::cerr << " Error in readMpp: File " << fileName
164  << " not found or locked" << std::endl;
165  std::abort();
166  }
167 
168  std::cout << "Reading mesh++ file" << std::endl;
169 
170  if ( ! readMppFileHead ( hstream, numberVertices, numberBoundaryVertices,
171  numberBoundaryFaces, numberBoundaryEdges, numberVolumes ) )
172  {
173  std::cerr << " Error While reading mesh++ file headers" << std::endl;
174  std::abort() ;
175  }
176 
177  hstream.close();
178 
179  UInt numberStoredEdges = numberBoundaryEdges;
180 
181  //Reopen the stream: I know it is stupid but this is how it goes
182  std::ifstream myStream ( fileName.c_str() );
183 
184  if ( myStream.fail() )
185  {
186  std::cerr << " Error in readMpp: file " << fileName
187  << " not found or locked" << std::endl;
188  std::abort();
189  }
190 
191  // Euler formulas
192  numberFaces = 2 * numberVolumes + ( numberBoundaryFaces / 2 );
193  numberEdges = numberVolumes + numberVertices +
194  ( 3 * numberBoundaryFaces - 2 * numberBoundaryVertices ) / 4;
195 
196  // Be a little verbose
197  if ( GeoShape::S_numPoints > 4 )
198  {
199  std::cout << "Quadratic Tetra Mesh (from Linear geometry)"
200  << std::endl;
201  numberPoints = numberVertices + numberEdges;
202  numberBoundaryPoints = numberBoundaryVertices + numberBoundaryEdges;
203  }
204 
205  else
206  {
207  std::cout << "Linear Tetra Mesh" << std::endl;
208  numberPoints = numberVertices;
209  numberBoundaryPoints = numberBoundaryVertices;
210  }
211 
212  std::cout << "Number of Vertices = " << std::setw ( 10 ) << numberVertices << std::endl
213  << "Number of Boundary Vertices = " << std::setw ( 10 ) << numberBoundaryVertices << std::endl;
214  oStr << "Number of Faces = " << std::setw ( 10 ) << numberFaces << std::endl
215  << "Number of Boundary Faces = " << std::setw ( 10 ) << numberBoundaryFaces << std::endl
216  << "Number of Edges = " << std::setw ( 10 ) << numberEdges << std::endl
217  << "Number of Boundary Edges = " << std::setw ( 10 ) << numberBoundaryEdges << std::endl;
218  std::cout << "Number of Points = " << std::setw ( 10 ) << numberPoints << std::endl
219  << "Number of Boundary Points = " << std::setw ( 10 ) << numberBoundaryPoints << std::endl
220  << "Number of Volumes = " << std::setw ( 10 ) << numberVolumes << std::endl;
221 
222  // Set all basic data structure:
223 
224  // I store all Points
225  mesh.setMaxNumPoints ( numberPoints, true );
226  mesh.setMaxNumGlobalPoints ( numberPoints );
227  mesh.setNumBPoints ( numberBoundaryPoints );
228  mesh.setNumVertices ( numberVertices );
229  mesh.setNumGlobalVertices ( numberVertices );
230  mesh.setNumBVertices ( numberBoundaryVertices );
231 
232  // Only Boundary Edges
233  mesh.setMaxNumEdges ( numberEdges );
234  mesh.setMaxNumGlobalEdges ( numberEdges );
235  mesh.setNumEdges ( numberEdges ); // Here the REAL number of edges (all of them)
236  mesh.setNumBEdges ( numberBoundaryEdges );
237 
238  // Only Boundary Faces
239  mesh.setMaxNumFaces ( numberBoundaryFaces );
240  mesh.setMaxNumGlobalFaces ( numberBoundaryFaces );
241  mesh.setNumFaces ( numberFaces ); // Here the REAL number of edges (all of them)
242  mesh.setNumBFaces ( numberBoundaryFaces );
243 
244  mesh.setMaxNumVolumes ( numberVolumes, true );
245  mesh.setMaxNumGlobalVolumes ( numberVolumes );
246 
247  mesh.setMarkerID ( regionFlag ); // Mark the region
248 
249 
250  typename mesh_Type::point_Type* pointerPoint = 0;
251  typename mesh_Type::edge_Type* pointerEdge = 0;
252  typename mesh_Type::face_Type* pointerFace = 0;
253  typename mesh_Type::volume_Type* pointerVolume = 0;
254 
255  // addPoint(), Face() and Edge() return a reference to the last stored point
256  // I use that information to set all point info, by using a pointer.
257 
258  UInt count = 0;
259  Int ibc;
260 
261  while ( nextGoodLine ( myStream, line ).good() )
262  {
263  if ( line.find ( "odes" ) != std::string::npos )
264  {
265 
266  std::string node_s = line.substr ( line.find_last_of ( ":" ) + 1 );
267  // _numberVertices=atoi(node_s);
268 
269  for ( i = 0; i < numberVertices; i++ )
270  {
271 #ifdef OLDMPPFILE
272  myStream >> x >> y >> z >> ity >> ibc;
273 #else
274 
275  myStream >> x >> y >> z >> ity >> ity_id;
276  if ( ity != 3 )
277  {
278  myStream >> ibc;
279  }
280 #endif
281 
282  if ( ity != 3 )
283  {
284  pointerPoint = &mesh.addPoint ( true, true );
285  pointerPoint->setId ( count );
286  ++count;
287 
288  //Boundary point. Boundary switch set by the mesh method.
289 
290  pointerPoint->setMarkerID ( markerID_Type ( ibc ) );
291  }
292  else
293  {
294  pointerPoint = &mesh.addPoint ( false, true );
295  pointerPoint->setId ( mesh.pointList.size() - 1 );
296  }
297 
298  pointerPoint->x() = x;
299  pointerPoint->y() = y;
300  pointerPoint->z() = z;
301  pointerPoint->setMarkerID ( markerID_Type ( ibc ) );
302  }
303 
304  oStr << "Vertices Read " << std::endl;
305  done++;
306 
307  if ( count != numberBoundaryVertices )
308  {
309  std::cerr << "NumB points inconsistent!" << std::endl;
310  }
311  }
312  if ( line.find ( "iangular" ) != std::string::npos )
313  {
314  oStr << "Reading Bfaces " << std::endl;
315  std::string node_s = line.substr ( line.find_last_of ( ":" ) + 1 );
316 
317  for ( i = 0; i < numberBoundaryFaces; i++ )
318  {
319 
320 #ifdef OLDMPPFILE
321  myStream >> p1 >> p2 >> p3 >> ity >> ibc;
322 #else
323 
324  myStream >> p1 >> p2 >> p3 >> ity >> ity_id >> ibc;
325 #endif
326  p1 -= idOffset;
327  p2 -= idOffset;
328  p3 -= idOffset; //get the 0-based numbering
329 
330  pointerFace = & ( mesh.addFace ( true ) ); // Only boundary faces
331 
332  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
333  pointerFace->setId ( i );
334  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
335  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
336  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
337  }
338 
339  oStr << "Boundary faces read " << std::endl;
340  done++;
341  }
342 
343  if ( line.find ( "Sides" ) != std::string::npos )
344  {
345  oStr << "Reading boundary edges " << std::endl;
346  std::string node_s = line.substr ( line.find_last_of ( ":" ) + 1 );
347  //_numberBoundaryEdges=atoi(node_s);
348 
349  for ( i = 0; i < numberStoredEdges; i++ )
350  {
351 #ifdef OLDMPPFILE
352  myStream >> p1 >> p2 >> ity >> ibc;
353 #else
354  myStream >> p1 >> p2 >> ity >> ity_id >> ibc;
355 #endif
356  p1 -= idOffset;
357  p2 -= idOffset; //get the 0-based numbering
358 
359  pointerEdge = &mesh.addEdge ( true ); // Only boundary edges.
360  pointerEdge->setMarkerID ( markerID_Type ( ibc ) );
361  pointerEdge->setId ( i );
362  pointerEdge->setPoint ( 0, mesh.point ( p1 ) ); // set edge conn.
363  pointerEdge->setPoint ( 1, mesh.point ( p2 ) ); // set edge conn.
364  }
365 
366  oStr << "Boundary edges read " << std::endl;
367  done++;
368  }
369 
370  count = 0;
371 
372  if ( line.find ( "etrahedral" ) != std::string::npos )
373  {
374  oStr << "Reading volumes " << std::endl;
375  std::string node_s = line.substr ( line.find_last_of ( ":" ) + 1 );
376 
377  for ( i = 0; i < numberVolumes; i++ )
378  {
379  myStream >> p1 >> p2 >> p3 >> p4;
380  p1 -= idOffset;
381  p2 -= idOffset;
382  p3 -= idOffset;
383  p4 -= idOffset; //get the 0-based numbering
384  pointerVolume = &mesh.addVolume();
385  pointerVolume->setId ( i );
386  pointerVolume->setPoint ( 0, mesh.point ( p1 ) );
387  pointerVolume->setPoint ( 1, mesh.point ( p2 ) );
388  pointerVolume->setPoint ( 2, mesh.point ( p3 ) );
389  pointerVolume->setPoint ( 3, mesh.point ( p4 ) );
390  count++;
391  }
392 
393  oStr << count << " Volume elements read" << std::endl;
394  done++;
395  }
396  }
397  // This part is to build a P2 mesh from a P1 geometry
398 
399  if ( GeoShape::S_numPoints > 4 )
400  {
401  MeshUtility::p2MeshFromP1Data ( mesh );
402  }
403 
404  myStream.close();
405 
406  // Test mesh
407  Switch sw;
408 
409  ///// CORRECTION JFG
410  //if (mesh.check(1, true,true))done=0;
411 
412  if ( !checkMesh3D ( mesh, sw, true, verbose, oStr, std::cerr, oStr ) )
413  {
414  std::abort(); // CORRECTION JFG
415  }
416 
417  Real vols[ 3 ];
418  getVolumeFromFaces ( mesh, vols, oStr );
419 
420  oStr << "Volume enclosed by the mesh computed by integration on boundary faces " << std::endl;
421  oStr << "INT(X) INT(Y) INT(Z) <- they should be equal and equal to" << std::endl
422  << " the voulume enclosed by the mesh " << std::endl;
423  oStr << vols[ 0 ] << " " << vols[ 1 ] << " " << vols[ 2 ] << std::endl;
424  oStr << "Boundary faces are defining a closed surface if "
425  << testClosedDomain ( mesh, oStr ) << std::endl
426  << "is (almost) zero" << std::endl;
427 
428  return done == 4 ;
429 }// Function readMppFile
430 
431 // ===================================================
432 // INRIA mesh readers
433 // ===================================================
434 
435 //! nextIntINRIAMeshField -
436 /*!
437  It gets an integer field from the std::string line
438  if it is not empty, otherwise from the input stream.
439  It assumes that the std::string is either empty or
440  it contains and integer. No check is made to verify this.
441 
442  @param line, the mesh data structure to fill in.
443  @param myStream, the name of the mesh file to read.
444  @return true if everything went fine, false otherwise.
445 */
446 
447 Int
448 nextIntINRIAMeshField ( std::string const& line,
449  std::istream& myStream );
450 
451 //! readINRIAMeshFileHead - It Reads all basic info from INRIA MESH.
452 /*!
453  It Reads all basic info from INRIA MESH file
454  so as to be able to properly dimension all arrays
455 
456  @param myStream,
457  @param numberVertices,
458  @param numberBoundaryVertices,
459  @param numberBoundaryFaces,
460  @param numberBoundaryEdges,
461  @param numberVolumes,
462  @param numberStoredFaces,
463  @param shape,
464  @param iSelect,
465  @return false if mesh check is unsuccessfull.
466 */
467 
468 bool
469 readINRIAMeshFileHead ( std::ifstream& myStream,
470  UInt& numberVertices,
471  UInt& numberBoundaryVertices,
472  UInt& numberBoundaryFaces,
473  UInt& numberBoundaryEdges,
474  UInt& numberVolumes,
475  UInt& numberStoredFaces,
476  ReferenceShapes& shape,
477  InternalEntitySelector iSelect = InternalEntitySelector() );
478 
479 //! readINRIAMeshFile - reads mesh++ Tetra meshes.
480 /*!
481  It converts Tetra meshes into Quadratic Tetra if needed it.
482 
483  @param mesh, the mesh data structure to fill in.
484  @param fileName, the name of the mesh file to read.
485  @param regionFlag, the identifier for the region.
486  @param verbose, setting it as true, the output is verbose (the default is false).
487  @param iSelect,
488  @return true if everything went fine, false otherwise.
489 */
490 
491 template <typename GeoShape, typename MC>
492 bool
493 readINRIAMeshFile ( RegionMesh<GeoShape, MC>& mesh,
494  std::string const& fileName,
495  markerID_Type regionFlag,
496  bool verbose = false,
497  InternalEntitySelector iSelect = InternalEntitySelector() )
498 {
499  typedef RegionMesh<GeoShape, MC> mesh_Type;
500 
501  const int idOffset = 1; //IDs in INRIA meshes start from 1
502 
503  std::string line;
504 
505 
506  Real x, y, z;
507 
508  UInt done = 0;
509  UInt i;
510  UInt numberVertices ( 0 ), numberBoundaryVertices ( 0 ),
511  numberFaces ( 0 ), numberBoundaryFaces ( 0 ),
512  numberPoints ( 0 ), numberBoundaryPoints ( 0 ),
513  numberEdges ( 0 ), numberBoundaryEdges ( 0 );
514  UInt numberVolumes ( 0 );
515  UInt numberStoredFaces ( 0 );
516  UInt p1, p2, p3, p4, p5, p6, p7, p8;
517 
518  std::stringstream discardedLog;
519 
520  std::vector<FiveNumbers> faceHelp;
521 
522  ReferenceShapes shape ( NONE );
523 
524  typename std::vector<FiveNumbers>::iterator faceHelpIterator;
525 
526  std::ostream& oStr = verbose ? std::cout : discardedLog;
527 
528  // open stream to read header
529 
530  std::ifstream hstream ( fileName.c_str() );
531 
532  if ( verbose )
533  {
534  std::cout << "Reading form file " << fileName << std::endl;
535  }
536 
537  if ( hstream.fail() )
538  {
539  std::cerr << " Error in readINRIAMeshFile = file " << fileName
540  << " not found or locked" << std::endl;
541  std::abort();
542  }
543 
544  if ( verbose )
545  {
546  std::cout << "Reading INRIA mesh file" << fileName << std::endl;
547  }
548 
549  if ( ! readINRIAMeshFileHead ( hstream, numberVertices, numberBoundaryVertices,
550  numberBoundaryFaces, numberBoundaryEdges,
551  numberVolumes, numberStoredFaces, shape, iSelect) )
552  {
553  std::cerr << " Error while reading INRIA mesh file headers" << std::endl;
554  std::abort() ;
555  }
556  //! Fix in case mesh file contains only a subset of the edges
557  UInt numberStoredEdges (numberBoundaryEdges);
558  hstream.close();
559 
560  //Reopen the stream: I know it is stupid but this is how it goes
561  std::ifstream myStream ( fileName.c_str() );
562 
563  if ( myStream.fail() )
564  {
565  std::cerr << " Error in readINRIAMeshFile = file " << fileName
566  << " not found or locked" << std::endl;
567  std::abort();
568  }
569 
570  ASSERT_PRE0 ( GeoShape::S_shape == shape, "INRIA Mesh file and mesh element shape is not consistent" );
571 
572  // Euler formulas to get number of faces and number of edges
573  numberFaces = 2 * numberVolumes + ( numberBoundaryFaces / 2 );
574  Int num1 = numberVertices + numberVolumes;
575  Int num2 = numberBoundaryVertices;
576  Int num3 = numberBoundaryFaces;
577 
578  numberEdges = ( 3 * num3 - 2 * num2 ) / 4 + num1;
579 
580  // numberEdges = (int) numberVolumes + numberVertices + ( 3 * numberBoundaryFaces + dummy ) / 4;
581 
582  // Be a little verbose
583  switch ( shape )
584  {
585 
586  case HEXA:
587  ASSERT_PRE0 ( GeoShape::S_numPoints == 8, "Sorry I can read only bilinear Hexa meshes" );
588  std::cout << "Linear Hexa mesh" << std::endl;
589  numberPoints = numberVertices;
590  numberBoundaryPoints = numberBoundaryVertices;
591  break;
592 
593  case TETRA:
594  if ( GeoShape::S_numPoints > 4 )
595  {
596  //if (GeoShape::S_numPoints ==6 )
597  std::cout << "Quadratic Tetra mesh (from linear geometry)" << std::endl;
598  numberPoints = numberVertices + numberEdges;
599 
600  /*
601  numberBoundaryPoints=numberBoundaryVertices+numberBoundaryEdges;
602  FALSE : numberBoundaryEdges is not known at this stage in a INRIA file
603  (JFG 07/2002)
604  I use the relation numberBoundaryVertices + numberBoundaryFaces - 2 = numberBoundaryEdges,
605  But, is it general (hole...)
606  (JFG 07/2002)
607 
608  numberBoundaryEdges = ( Int ( numberBoundaryVertices + numberBoundaryFaces - Int (2 ) )
609  > 0 ? ( numberBoundaryVertices + numberBoundaryFaces - 2 ) : 0 );
610  numberBoundaryPoints = ( Int ( numberBoundaryVertices +
611  ( numberBoundaryVertices + numberBoundaryFaces - Int ( 2 ) ) )
612  > 0 ?numberBoundaryVertices +
613  ( numberBoundaryVertices + numberBoundaryFaces - 2 ) : 0 );
614  */
615  }
616 
617  else
618  {
619  if ( verbose )
620  {
621  std::cout << "Linear Tetra Mesh" << std::endl;
622  }
623 
624  numberPoints = numberVertices;
625  numberBoundaryPoints = numberBoundaryVertices;
626  numberBoundaryEdges = ( Int ( numberBoundaryVertices + numberBoundaryFaces - Int ( 2 ) )
627  > 0 ? ( numberBoundaryVertices + numberBoundaryFaces - 2 ) : 0 );
628  }
629 
630  break;
631 
632  default:
633  ERROR_MSG ( "Current version of INRIA Mesh file reader only accepts TETRA and HEXA" );
634  }
635 
636  oStr << "Number of Vertices = " << std::setw ( 10 ) << numberVertices << std::endl
637  << "Number of BVertices = " << std::setw ( 10 ) << numberBoundaryVertices << std::endl
638  << "Number of Faces = " << std::setw ( 10 ) << numberFaces << std::endl
639  << "Number of Boundary Faces = " << std::setw ( 10 ) << numberBoundaryFaces << std::endl
640  << "Number of Stored Faces = " << std::setw ( 10 ) << numberStoredFaces << std::endl
641  << "Number of Edges = " << std::setw ( 10 ) << numberEdges << std::endl
642  << "Number of Boundary Edges = " << std::setw ( 10 ) << numberBoundaryEdges << std::endl
643  << "Number of Stored Edges = " << std::setw ( 10 ) << numberStoredEdges << std::endl
644  << "Number of Points = " << std::setw ( 10 ) << numberPoints << std::endl
645  << "Number of Boundary Points = " << std::setw ( 10 ) << numberBoundaryPoints << std::endl
646  << "Number of Volumes = " << std::setw ( 10 ) << numberVolumes << std::endl;
647 
648  // Set all basic data structure
649 
650  // I store all Points
651  mesh.setMaxNumPoints ( numberPoints, true );
652  mesh.setMaxNumGlobalPoints ( numberPoints );
653  mesh.setNumBPoints ( numberBoundaryPoints );
654  mesh.setNumVertices ( numberVertices );
655  mesh.setNumGlobalVertices ( numberVertices );
656  mesh.setNumBVertices ( numberBoundaryVertices );
657  // Only Boundary Edges (in a next version I will allow for different choices)
658  mesh.setMaxNumEdges ( numberBoundaryEdges );
659  mesh.setMaxNumGlobalEdges ( numberEdges );
660  mesh.setNumEdges ( numberEdges ); // Here the REAL number of edges (all of them)
661  mesh.setNumBEdges ( numberBoundaryEdges );
662  // Only Boundary Faces
663  mesh.setMaxNumFaces ( numberStoredFaces );
664  mesh.setMaxNumGlobalFaces ( numberBoundaryFaces );
665  mesh.setNumFaces ( numberFaces ); // Here the REAL number of faces (all of them)
666  mesh.setNumBFaces ( numberBoundaryFaces );
667 
668  mesh.setMaxNumVolumes ( numberVolumes, true );
669  mesh.setMaxNumGlobalVolumes ( numberVolumes );
670 
671  mesh.setMarkerID ( regionFlag ); // Add Marker to list of Markers
672 
673  typedef typename mesh_Type::point_Type point_Type;
674  typedef typename mesh_Type::volume_Type volume_Type;
675 
676 
677  point_Type* pointerPoint = 0;
678  typename mesh_Type::edge_Type* pointerEdge = 0;
679  typename mesh_Type::face_Type* pointerFace = 0;
680  volume_Type* pointerVolume = 0;
681  // addPoint()/Face()/Edge() returns a reference to the last stored point
682  // I use that information to set all point info, by using a pointer.
683 
684  UInt count = 0;
685  Int ibc;
686 
687  // To account for internal faces
688  if ( numberStoredFaces > numberBoundaryFaces )
689  {
690  faceHelp.resize ( numberStoredFaces - numberBoundaryFaces );
691  faceHelpIterator = faceHelp.begin();
692 
693  oStr << "WARNING: The mesh file (apparently) contains "
694  << numberStoredFaces - numberBoundaryFaces << " internal faces" << std::endl;
695 
696  }
697 
698  while ( nextGoodLine ( myStream, line ).good() )
699  {
700  if ( line.find ( "Vertices" ) != std::string::npos )
701  {
702  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "s" ) + 1 ), myStream );
703 
704  for ( i = 0; i < numberVertices; i++ )
705  {
706  myStream >> x >> y >> z >> ibc;
707 
708  if ( !iSelect (markerID_Type (ibc) ) )
709  {
710  ++count;
711  // Boundary point. Boundary switch set by the mesh method.
712  pointerPoint = &mesh.addPoint ( true, true );
713  pointerPoint->setMarkerID ( markerID_Type ( ibc ) );
714  }
715 
716  else
717  {
718  pointerPoint = &mesh.addPoint ( false, true );
719  }
720 
721  pointerPoint->setId ( i );
722  pointerPoint->x() = x;
723  pointerPoint->y() = y;
724  pointerPoint->z() = z;
725  pointerPoint->setMarkerID ( markerID_Type ( ibc ) );
726  }
727 
728  oStr << "Vertices read " << std::endl;
729  oStr << "Size of the node storage is "
730  << count* sizeof ( point_Type ) / 1024. / 1024. << std::endl;
731  done++;
732 
733  if ( count != numberBoundaryVertices )
734  {
735  std::cerr << "Number boundary points inconsistent!" << std::endl;
736  }
737  }
738 
739  if ( line.find ( "Triangles" ) != std::string::npos )
740  {
741  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "s" ) + 1 ), myStream );
742  oStr << "Reading boundary faces " << std::endl;
743 
744  for ( i = 0; i < numberStoredFaces; i++ )
745  {
746  myStream >> p1 >> p2 >> p3 >> ibc;
747  p1 -= idOffset;
748  p2 -= idOffset;
749  p3 -= idOffset; //get the 0-based numbering
750 
751  if ( numberStoredFaces > numberBoundaryFaces )
752  {
753  if ( mesh.point ( p1 ).boundary() && mesh.point ( p2 ).boundary() &&
754  mesh.point ( p3 ).boundary() )
755  {
756  pointerFace = & ( mesh.addFace ( true ) ); // Boundary faces
757  pointerFace->setId ( i );
758  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
759  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
760  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
761  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
762 
763  }
764 
765  else
766  {
767  faceHelpIterator->i1 = p1;
768  faceHelpIterator->i2 = p2;
769  faceHelpIterator->i3 = p3;
770  faceHelpIterator->ibc = ibc;
771 
772  ++faceHelpIterator;
773  }
774  }
775 
776  else
777  {
778 
779  pointerFace = & ( mesh.addFace ( true ) ); // Only boundary faces
780 
781  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
782  pointerFace->setId ( i );
783  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
784  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
785  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
786  }
787  }
788 
789  UInt extraFaceCounter = numberBoundaryFaces;
790  for ( faceHelpIterator = faceHelp.begin();
791  faceHelpIterator != faceHelp.end(); ++faceHelpIterator, extraFaceCounter++ )
792  {
793  p1 = faceHelpIterator->i1;
794  p2 = faceHelpIterator->i2;
795  p3 = faceHelpIterator->i3;
796  ibc = faceHelpIterator->ibc;
797  pointerFace = & ( mesh.addFace ( false ) ); // INTERNAL FACE
798  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
799  pointerFace->setId ( extraFaceCounter );
800  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
801  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
802  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
803  }
804 
805  oStr << "Boundary faces read " << std::endl;
806  done++;
807  }
808 
809  if ( line.find ( "Quadrilaterals" ) != std::string::npos )
810  {
811  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "s" ) + 1 ), myStream );
812 
813  oStr << "Reading boundary faces " << std::endl;
814 
815  for (UInt i = 0; i < numberBoundaryFaces; i++ )
816  {
817  myStream >> p1 >> p2 >> p3 >> p4 >> ibc;
818  p1 -= idOffset;
819  p2 -= idOffset;
820  p3 -= idOffset;
821  p4 -= idOffset; //get the 0-based numbering
822 
823  if ( numberStoredFaces > numberBoundaryFaces )
824  {
825  if ( mesh.point ( p1 ).boundary() && mesh.point ( p2 ).boundary() &&
826  mesh.point ( p3 ).boundary() )
827  {
828  pointerFace = & ( mesh.addFace ( true ) ); // Boundary faces
829  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
830  pointerFace->setId ( i );
831  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
832  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
833  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
834  pointerFace->setPoint ( 3, mesh.point ( p4 ) ); // set face conn.
835 
836  }
837 
838  else
839  {
840  faceHelpIterator->i1 = p1;
841  faceHelpIterator->i2 = p2;
842  faceHelpIterator->i3 = p3;
843  faceHelpIterator->i4 = p4;
844  faceHelpIterator->ibc = ibc;
845 
846  ++faceHelpIterator;
847  }
848  }
849 
850  else
851  {
852  pointerFace = & ( mesh.addFace ( true ) ); // Only boundary faces
853  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
854  pointerFace->setId ( i );
855  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
856  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
857  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
858  pointerFace->setPoint ( 3, mesh.point ( p4 ) ); // set face conn.
859  }
860  }
861 
862  oStr << "Boundary faces read " << std::endl;
863 
864  UInt extraFaceCounter = numberBoundaryFaces;
865  for ( faceHelpIterator = faceHelp.begin();
866  faceHelpIterator != faceHelp.end(); ++faceHelpIterator, extraFaceCounter++ )
867  {
868  p1 = faceHelpIterator->i1;
869  p2 = faceHelpIterator->i2;
870  p3 = faceHelpIterator->i3;
871  p4 = faceHelpIterator->i4;
872  ibc = faceHelpIterator->ibc;
873  pointerFace = & ( mesh.addFace ( false ) ); // INTERNAL FACE
874  pointerFace->setMarkerID ( markerID_Type ( ibc ) );
875  pointerFace->setId ( extraFaceCounter );
876  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
877  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
878  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
879  pointerFace->setPoint ( 3, mesh.point ( p4 ) ); // set face conn.
880  }
881  done++;
882  }
883 
884  if ( line.find ( "Edges" ) != std::string::npos )
885  {
886  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "s" ) + 1 ), myStream );
887  oStr << "Reading stored edges " << std::endl;
888 
889  for ( i = 0; i < numberStoredEdges; i++ )
890  {
891  myStream >> p1 >> p2 >> ibc;
892  p1 -= idOffset;
893  p2 -= idOffset; //get the 0-based numbering
894  pointerEdge = &mesh.addEdge ( true ); // Only boundary edges.
895  pointerEdge->setMarkerID ( markerID_Type ( ibc ) );
896  pointerEdge->setId ( i );
897  pointerEdge->setPoint ( 0, mesh.point ( p1 ) ); // set edge conn.
898  pointerEdge->setPoint ( 1, mesh.point ( p2 ) ); // set edge conn.
899  }
900  oStr << " Stored edges read " << std::endl;
901  done++;
902  }
903 
904  if ( line.find ( "Tetrahedra" ) != std::string::npos )
905  {
906  count = 0;
907  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "a" ) + 1 ), myStream );
908  oStr << "Reading volumes " << std::endl;
909 
910  for ( i = 0; i < numberVolumes; i++ )
911  {
912  myStream >> p1 >> p2 >> p3 >> p4 >> ibc;
913  p1 -= idOffset;
914  p2 -= idOffset;
915  p3 -= idOffset;
916  p4 -= idOffset; //get the 0-based numbering
917  pointerVolume = &mesh.addVolume();
918  pointerVolume->setId ( i );
919  pointerVolume->setPoint ( 0, mesh.point ( p1 ) );
920  pointerVolume->setPoint ( 1, mesh.point ( p2 ) );
921  pointerVolume->setPoint ( 2, mesh.point ( p3 ) );
922  pointerVolume->setPoint ( 3, mesh.point ( p4 ) );
923  pointerVolume->setMarkerID ( markerID_Type ( ibc ) );
924  count++;
925  }
926  oStr << "size of the volume storage is " << sizeof ( volume_Type ) * count / 1024. / 1024.
927  << " Mo." << std::endl;
928  oStr << count << " Volume elements read" << std::endl;
929  done++;
930  }
931 
932  if ( line.find ( "Hexahedra" ) != std::string::npos )
933  {
934  count = 0;
935  nextIntINRIAMeshField ( line.substr ( line.find_last_of ( "a" ) + 1 ), myStream );
936  oStr << "Reading volumes " << std::endl;
937  for ( i = 0; i < numberVolumes; i++ )
938  {
939  myStream >> p1 >> p2 >> p3 >> p4 >> p5 >> p6 >> p7 >> p8 >> ibc;
940  p1 -= idOffset;
941  p2 -= idOffset;
942  p3 -= idOffset;
943  p4 -= idOffset; //get the 0-based numbering
944  p5 -= idOffset;
945  p6 -= idOffset;
946  p7 -= idOffset;
947  p8 -= idOffset; //get the 0-based numbering
948  pointerVolume = &mesh.addVolume();
949  pointerVolume->setId ( i );
950  pointerVolume->setPoint ( 0, mesh.point ( p1 ) );
951  pointerVolume->setPoint ( 1, mesh.point ( p2 ) );
952  pointerVolume->setPoint ( 2, mesh.point ( p3 ) );
953  pointerVolume->setPoint ( 3, mesh.point ( p4 ) );
954  pointerVolume->setPoint ( 4, mesh.point ( p5 ) );
955  pointerVolume->setPoint ( 5, mesh.point ( p6 ) );
956  pointerVolume->setPoint ( 6, mesh.point ( p7 ) );
957  pointerVolume->setPoint ( 7, mesh.point ( p8 ) );
958  pointerVolume->setMarkerID ( markerID_Type ( ibc ) );
959 
960  count++;
961  }
962  oStr << count << " Volume elements read" << std::endl;
963  done++;
964  }
965 
966  }
967 
968  // Test the mesh
969  Switch sw;
970 
971  // this if is the verbose version
972  if ( !checkMesh3D ( mesh, sw, true, verbose, oStr, std::cerr, oStr ) )
973  {
974  std::abort();
975  }
976 
977  // This part is to build a P2 mesh from a P1 geometry
978  if ( shape == TETRA && GeoShape::S_numPoints > 4 )
979  {
980  MeshUtility::p2MeshFromP1Data ( mesh );
981  }
982 
983  myStream.close();
984 
985  Real vols[ 3 ];
986  getVolumeFromFaces ( mesh, vols, oStr );
987 
988  oStr << "Volume enclosed by the mesh computed by integration on boundary faces" << std::endl;
989  oStr << "INT(X) INT(Y) INT(Z) <- they should be equal and equal to" << std::endl
990  << " the voulume enclosed by the mesh" << std::endl;
991  oStr << vols[ 0 ] << " " << vols[ 1 ] << " " << vols[ 2 ] << std::endl;
992  oStr << "Boundary faces are defining a closed surface if "
993  << testClosedDomain ( mesh, oStr ) << std::endl
994  << "is (almost) zero" << std::endl;
995 
996  return done == 4 ;
997 }// Function readINRIAMeshFile
998 
999 // ===================================================
1000 // GMSH mesh readers
1001 // ===================================================
1002 
1003 //! readGmshFile - it reads a GMSH mesh file
1004 /*!
1005  It reads a 3D gmsh mesh file and store it in a RegionMesh.
1006 
1007  @param mesh mesh data structure to fill in
1008  @param fileName name of the gmsh mesh file to read
1009  @param regionFlag identifier for the region
1010  @param verbose whether the function shall be verbose
1011  @return true if everything went fine, false otherwise
1012 */
1013 
1014 template <typename GeoShape, typename MC>
1015 bool
1016 readGmshFile ( RegionMesh<GeoShape, MC>& mesh,
1017  const std::string& fileName,
1018  markerID_Type regionFlag,
1019  bool verbose = false )
1020 {
1021  typedef RegionMesh<GeoShape, MC> mesh_Type;
1022 
1023  const int idOffset = 1; //IDs in GMESH files start from 1
1024 
1025  std::ifstream inputFile ( fileName.c_str() );
1026 
1027 #ifdef HAVE_LIFEV_DEBUG
1028  debugStream ( 8000 ) << "Gmsh reading: " << fileName << "\n";
1029 #endif
1030 
1031  // char buffer[256];
1032  std::string buffer;
1033 
1034  for (Int ii = 0; ii < 6; ++ii)
1035  {
1036  inputFile >> buffer;
1037  std::cout << "buffer = " << buffer << "\n";
1038  }
1039 
1040  UInt numberNodes;
1041  inputFile >> numberNodes;
1042 
1043 #ifdef HAVE_LIFEV_DEBUG
1044  debugStream ( 8000 ) << "Number of nodes = " << numberNodes;
1045 #endif
1046 
1047  // Add Marker to list of Markers
1048  mesh.setMarkerID ( regionFlag );
1049 
1050 
1051  std::vector<Real> x ( 3 * numberNodes );
1052  std::vector<bool> isonboundary ( numberNodes );
1053  std::vector<UInt> whichboundary ( numberNodes );
1054 
1055 #ifdef HAVE_LIFEV_DEBUG
1056  debugStream ( 8000 ) << "Reading " << numberNodes << " nodes\n";
1057 #endif
1058 
1059  std::map<Int, Int> itoii;
1060 
1061  for ( UInt i = 0; i < numberNodes; ++i )
1062  {
1063  UInt ni;
1064  inputFile >> ni
1065  >> x[ 3 * i ]
1066  >> x[ 3 * i + 1 ]
1067  >> x[ 3 * i + 2 ];
1068 
1069  itoii[ ni - idOffset] = i;
1070  }
1071  inputFile >> buffer;
1072 
1073 #ifdef HAVE_LIFEV_DEBUG
1074  debugStream ( 8000 ) << "buffer = " << buffer << "\n";
1075 #endif
1076 
1077  inputFile >> buffer;
1078 
1079 #ifdef HAVE_LIFEV_DEBUG
1080  debugStream ( 8000 ) << "buffer = " << buffer << "\n";
1081 #endif
1082 
1083  UInt numberElements;
1084  inputFile >> numberElements;
1085 
1086  typename mesh_Type::edge_Type* pointerEdge = 0;
1087  typename mesh_Type::face_Type* pointerFace = 0;
1088  typename mesh_Type::volume_Type* pointerVolume = 0;
1089 
1090 #ifdef HAVE_LIFEV_DEBUG
1091  debugStream ( 8000 ) << "number of elements: " << numberElements << "\n";
1092 #endif
1093 
1094  std::vector<std::vector<int> > e ( numberElements );
1095  std::vector<int> et ( numberElements );
1096  std::vector<int> etype ( numberElements );
1097  std::vector<int> gt ( 32 );
1098  gt.assign ( 32, 0 );
1099 
1100  for ( UInt i = 0; i < numberElements; ++i )
1101  {
1102  Int ne, t, np;
1103 
1104  inputFile >> buffer;
1105  inputFile >> ne;
1106 
1107  switch ( ne )
1108  {
1109  case (2) :
1110  np = 3;
1111  break;
1112 
1113  case (4) :
1114  np = 4;
1115  break;
1116 
1117  case (15) :
1118  np = 1;
1119  break;
1120 
1121  default:
1122  np = 0;
1123 
1124 #ifdef HAVE_LIFEV_DEBUG
1125  debugStream ( 8000 ) << "Element type unknown " << ne << "\n";
1126 #endif
1127 
1128  ASSERT ( true, "Elements type unsupported.\n" )
1129  }
1130 
1131  inputFile >> t;
1132 
1133  //std::debug() << t << " ";
1134 
1135  bool ibcSet = false;
1136  Int flag = 0;
1137  Int tag ( 0 );
1138 
1139  for ( Int iflag = 0; iflag < t; ++iflag )
1140  {
1141  inputFile >> flag;
1142 
1143  if ( !ibcSet )
1144  {
1145  tag = flag;
1146  ibcSet = true;
1147  }
1148  }
1149 
1150  ++gt[ ne ];
1151 
1152  etype[ i ] = ne;
1153 
1154  et[ i ] = tag;
1155  e[ i ].resize ( np );
1156 
1157  Int p = 0;
1158 
1159  while ( p != np )
1160  {
1161  Int node;
1162  inputFile >> node;
1163  e[ i ][ p ] = itoii[ node - idOffset ];
1164  ++p;
1165  }
1166  }
1167 
1168 
1169  // Euler formulas
1170  UInt n_volumes = gt[ 4 ];
1171  UInt n_faces_boundary = gt[ 2 ];
1172  UInt n_faces_total = 2 * n_volumes + ( n_faces_boundary / 2 );
1173 
1174  mesh.setMaxNumGlobalPoints ( numberNodes );
1175  // Only Boundary Edges (in a next version I will allow for different choices)
1176  mesh.setMaxNumEdges ( gt[ 1 ] );
1177  mesh.setNumEdges ( gt[ 1 ] ); // Here the REAL number of edges (all of them)
1178  mesh.setNumBEdges ( gt[ 1 ] );
1179  mesh.setMaxNumGlobalEdges ( gt[ 1 ] );
1180 
1181 #ifdef HAVE_LIFEV_DEBUG
1182  debugStream ( 8000 ) << "number of edges= " << gt[ 1 ] << "\n";
1183 #endif
1184 
1185  // Only Boundary Faces
1186  mesh.setMaxNumFaces ( n_faces_total );
1187  mesh.setNumFaces ( n_faces_total ); // Here the REAL number of edges (all of them)
1188  //mesh.setMaxNumFaces( n_faces_boundary );
1189  //mesh.setNumFaces ( n_faces_boundary ); // Here the REAL number of edges (all of them)
1190  mesh.setNumBFaces ( n_faces_boundary );
1191  mesh.setMaxNumGlobalFaces ( n_faces_total );
1192 
1193 #ifdef HAVE_LIFEV_DEBUG
1194  debugStream ( 8000 ) << "number of faces = " << n_faces_boundary << "\n";
1195 #endif
1196 
1197  mesh.setMaxNumVolumes ( n_volumes, true );
1198  mesh.setMaxNumGlobalVolumes ( n_volumes );
1199 
1200 #ifdef HAVE_LIFEV_DEBUG
1201  debugStream ( 8000 ) << "number of volumes = " << n_volumes << "\n";
1202 #endif
1203 
1204  isonboundary.assign ( numberNodes, false );
1205  whichboundary.assign ( numberNodes, 0 );
1206 
1207  for ( UInt i = 0; i < numberElements; ++i )
1208  {
1209  switch ( etype[ i ] )
1210  {
1211  // triangular faces (linear)
1212  case 2:
1213  {
1214  isonboundary[ e[ i ][ 0 ] ] = true;
1215  isonboundary[ e[ i ][ 1 ] ] = true;
1216  isonboundary[ e[ i ][ 2 ] ] = true;
1217 
1218  whichboundary[ e[ i ][ 0 ] ] = et[ i ];
1219  whichboundary[ e[ i ][ 1 ] ] = et[ i ];
1220  whichboundary[ e[ i ][ 2 ] ] = et[ i ];
1221  }
1222  }
1223  }
1224  // add the point to the mesh
1225  typename mesh_Type::point_Type* pointerPoint = 0;
1226 
1227  mesh.setMaxNumPoints ( numberNodes, true );
1228  mesh.setNumVertices ( numberNodes );
1229  mesh.setNumBVertices ( std::count ( isonboundary.begin(), isonboundary.end(), true ) );
1230  mesh.setNumBPoints ( mesh.numBVertices() );
1231 
1232 #ifdef HAVE_LIFEV_DEBUG
1233  debugStream ( 8000 ) << "number of points : " << mesh.numPoints() << "\n";
1234  debugStream ( 8000 ) << "number of boundary points : " << mesh.numBPoints() << "\n";
1235  debugStream ( 8000 ) << "number of vertices : " << mesh.numVertices() << "\n";
1236  debugStream ( 8000 ) << "number of boundary vertices : " << mesh.numBVertices() << "\n";
1237 #endif
1238 
1239  for ( UInt i = 0; i < numberNodes; ++i )
1240  {
1241  pointerPoint = &mesh.addPoint ( isonboundary[ i ], true );
1242  pointerPoint->setMarkerID ( whichboundary[ i ] );
1243  pointerPoint->setId ( i );
1244  pointerPoint->x() = x[ 3 * i ];
1245  pointerPoint->y() = x[ 3 * i + 1 ];
1246  pointerPoint->z() = x[ 3 * i + 2 ];
1247  }
1248 
1249  // add the element to the mesh
1250  for ( UInt i = 0; i < numberElements; ++i )
1251  {
1252  switch ( etype[ i ] )
1253  {
1254  // segment(linear)
1255  case 1:
1256  {
1257  pointerEdge = & ( mesh.addEdge ( true ) );
1258  pointerEdge->setMarkerID ( markerID_Type ( et[ i ] ) );
1259  pointerEdge->setId ( i );
1260  pointerEdge->setPoint ( 0, mesh.point ( e[ i ][ 0 ] ) );
1261  pointerEdge->setPoint ( 1, mesh.point ( e[ i ][ 1 ] ) );
1262 
1263 
1264 
1265  }
1266  break;
1267 
1268  // triangular faces (linear)
1269  case 2:
1270  {
1271  pointerFace = & ( mesh.addFace ( true ) );
1272  pointerFace->setMarkerID ( markerID_Type ( et[ i ] ) );
1273  pointerFace->setId ( i );
1274  pointerFace->setPoint ( 0, mesh.point ( e[ i ][ 0 ] ) );
1275  pointerFace->setPoint ( 1, mesh.point ( e[ i ][ 1 ] ) );
1276  pointerFace->setPoint ( 2, mesh.point ( e[ i ][ 2 ] ) );
1277 
1278  }
1279  break;
1280 
1281  // quadrangular faces(linear)
1282  case 3:
1283  {
1284  pointerFace = & ( mesh.addFace ( true ) );
1285  pointerFace->setMarkerID ( markerID_Type ( et[ i ] ) );
1286  pointerFace->setId ( i );
1287  pointerFace->setPoint ( 0, mesh.point ( e[ i ][ 0 ] ) );
1288  pointerFace->setPoint ( 1, mesh.point ( e[ i ][ 1 ] ) );
1289  pointerFace->setPoint ( 2, mesh.point ( e[ i ][ 2 ] ) );
1290  pointerFace->setPoint ( 3, mesh.point ( e[ i ][ 3 ] ) );
1291  }
1292  break;
1293 
1294  // tetrahedrons(linear)
1295  case 4:
1296  {
1297  pointerVolume = & ( mesh.addVolume() );
1298  pointerVolume->setId ( i );
1299  pointerVolume->setMarkerID ( markerID_Type ( et[ i ] ) );
1300  pointerVolume->setPoint ( 0, mesh.point ( e[ i ][ 0 ] ) );
1301  pointerVolume->setPoint ( 1, mesh.point ( e[ i ][ 1 ] ) );
1302  pointerVolume->setPoint ( 2, mesh.point ( e[ i ][ 2 ] ) );
1303  pointerVolume->setPoint ( 3, mesh.point ( e[ i ][ 3 ] ) );
1304  }
1305  break;
1306 
1307  // hexahedrons(linear)
1308  case 5:
1309  {
1310  pointerVolume = & ( mesh.addVolume() );
1311 
1312  pointerVolume->setId ( i );
1313  pointerVolume->setMarkerID ( markerID_Type ( et[ i ] ) );
1314  pointerVolume->setPoint ( 0, mesh.point ( e[ i ][ 0 ] ) );
1315  pointerVolume->setPoint ( 1, mesh.point ( e[ i ][ 1 ] ) );
1316  pointerVolume->setPoint ( 2, mesh.point ( e[ i ][ 2 ] ) );
1317  pointerVolume->setPoint ( 3, mesh.point ( e[ i ][ 3 ] ) );
1318  pointerVolume->setPoint ( 4, mesh.point ( e[ i ][ 4 ] ) );
1319  pointerVolume->setPoint ( 5, mesh.point ( e[ i ][ 5 ] ) );
1320  pointerVolume->setPoint ( 6, mesh.point ( e[ i ][ 6 ] ) );
1321  pointerVolume->setPoint ( 7, mesh.point ( e[ i ][ 7 ] ) );
1322  }
1323  break;
1324  }
1325  }
1326 
1327  Switch sw;
1328 
1329  std::stringstream discardedLog;
1330  std::ostream& oStr = verbose ? std::cout : discardedLog;
1331 
1332  if ( checkMesh3D (mesh, sw, true, verbose, oStr, std::cerr, oStr) == false )
1333  {
1334  std::ostringstream ex;
1335  ex << "invalid mesh from GSMH";
1336 
1337  throw std::logic_error ( ex.str() );
1338  }
1339 
1340  return true;
1341 } // Function readGmshFile
1342 
1343 // ===================================================
1344 // NetGen mesh readers
1345 // ===================================================
1346 
1347 //! readNetgenMesh - reads mesh++ Tetra meshes.
1348 /*!
1349  It reads a 3D NetGen mesh file and store it in a RegionMesh.
1350 
1351  @param mesh mesh data structure to fill in.
1352  @param fileName name of the gmsh mesh file to read.
1353  @param regionFlag identifier for the region.
1354  @param verbose whether the function shall be verbose.
1355  @return true if everything went fine, false otherwise.
1356 */
1357 
1358 template<typename GeoShape, typename MC>
1359 bool
1360 readNetgenMesh (RegionMesh<GeoShape, MC>& mesh,
1361  const std::string& fileName,
1362  markerID_Type regionFlag,
1363  bool verbose = false )
1364 {
1365  typedef RegionMesh<GeoShape, MC> mesh_Type;
1366 
1367  const int idOffset = 1; //IDs in netgen starts from 1
1368 
1369  // I will extract lines from iostream
1370  std::string line;
1371 
1372  // number of Geo Elements
1373  UInt numberVertices ( 0 ), numberBoundaryVertices ( 0 ),
1374  numberPoints ( 0 ), numberBoundaryPoints ( 0 ),
1375  numberEdges ( 0 ), numberBoundaryEdges ( 0 ),
1376  numberFaces ( 0 ), numberBoundaryFaces ( 0 ),
1377  numberVolumes ( 0 );
1378 
1379  // During the first access to file, build list of structures
1380  Vector pointCoordinates;
1381  std::vector<UInt> edgePointID, facePointID, volumePointID;
1382 
1383  // for poit i-th: is it a boundary point?
1384  std::vector<bool> boundaryPoint;
1385 
1386  // build a list of boundary edges, since netgen is not writing all of them
1387  MeshElementBareHandler<BareEdge> bihBedges;
1388 
1389  // flags for boundary entities
1390  std::vector<markerID_Type> bcnsurf, bcnpoints;
1391 
1392  // bitstream to check which file section has already been visited
1393  UInt flag;
1394 
1395  typename MC::pointMarker_Type pointMarker;
1396  typename MC::edgeMarker_Type edgeMarker;
1397  typename MC::faceMarker_Type faceMarker;
1398 
1399  // open file stream to look for points information
1400  std::ifstream fstreamp ( fileName.c_str() );
1401  if (fstreamp.fail() )
1402  {
1403  std::cerr << "Error in readNetgenMesh: File not found or locked" << std::endl;
1404  std::abort();
1405  }
1406 
1407  std::cout << "Reading netgen mesh file: " << fileName << std::endl;
1408  getline ( fstreamp, line );
1409 
1410  if ( line.find ( "mesh3d" ) == std::string::npos )
1411  {
1412  std::cerr << "Error in readNetgenMesh: mesh file is not in mesh3d format (netgen)"
1413  << std::endl;
1414  std::abort();
1415  }
1416 
1417  /*
1418  I assume as I tested that faces stored are only boundary faces
1419  and edges stored are only a part of boundary ones with
1420  inside a face with a common marker, instead points can be
1421  inside the domain too, but this format doesn't say me
1422  which, so I'll find them from myself
1423  */
1424 
1425  flag = 1 | 2 | 4 | 8;
1426 
1427  while ( getline ( fstreamp, line ) )
1428  {
1429 
1430  if ( line.find ( "points" ) != std::string::npos && flag & 1 )
1431  {
1432  getline ( fstreamp, line );
1433  std::stringstream parseLine ( line );
1434 
1435  parseLine >> numberVertices;
1436  std::cout << "[readNetgenMesh] found " << numberVertices << " vertices... " << std::flush;
1437 
1438  bcnpoints.resize ( numberVertices + 1 );
1439  boundaryPoint.resize ( numberVertices + 1 );
1440  pointCoordinates.resize ( nDimensions * numberVertices );
1441 
1442  for ( UInt i = 0; i < numberVertices; i++ )
1443  {
1444  // helping variables to read netgen file fields
1445  Real x, y, z;
1446 
1447  getline ( fstreamp, line );
1448  std::stringstream parseLine ( line );
1449 
1450  parseLine >> x >> y >> z;
1451 
1452  pointCoordinates[ i * nDimensions ] = x;
1453  pointCoordinates[ i * nDimensions + 1 ] = y;
1454  pointCoordinates[ i * nDimensions + 2 ] = z;
1455  boundaryPoint [ i ] = false;
1456  bcnpoints[ i ] = pointMarker.nullMarkerID();
1457  }
1458  boundaryPoint [ numberVertices ] = false;
1459  bcnpoints[ numberVertices ] = pointMarker.nullMarkerID();
1460 
1461  // done parsing point section
1462  flag &= ~1;
1463  std::cout << "loaded." << std::endl;
1464  break;
1465  }
1466  }
1467  fstreamp.close();
1468 
1469  /* as Forma told, I have found sometimes problems with seekg, becouse
1470  it seems that sometimes it put the stream in a non consistent state,
1471  so I open new files instead of seeking, I have had not compile errors
1472  only sometimes problems on big mesh files with seekg
1473  */
1474  // fstream.seekg(0,std::ios_base::beg);
1475 
1476  /* I assume as I tested that edges stored are only a part of
1477  boundary ones, so really they are meaningless to me
1478  unless I'm working with 2D meshes: in that case
1479  I can extract from here information about boundary vertices */
1480  std::ifstream fstreame ( fileName.c_str() );
1481  while ( getline ( fstreame, line ) )
1482  {
1483 
1484  if ( line.find ( "edgesegmentsgi2" ) != std::string::npos && flag & 2 )
1485  {
1486  getline ( fstreame, line );
1487  std::stringstream parseLine ( line );
1488 
1489  parseLine >> numberBoundaryEdges; // this will be discarded in 3D meshes
1490 
1491  std::cout << "[readNetgenMesh] found " << numberBoundaryEdges << " boundary edges... " << std::flush;
1492 
1493  edgePointID.resize ( 2 * numberBoundaryEdges );
1494 
1495  for ( UInt i = 0; i < numberBoundaryEdges; ++i)
1496  {
1497  UInt surfnr, t, p1, p2;
1498 
1499  getline ( fstreame, line );
1500  std::stringstream parseLine ( line );
1501 
1502  parseLine >> surfnr >> t >> p1 >> p2;
1503 
1504  edgePointID[ 2 * i ] = p1 - idOffset;
1505  edgePointID[ 2 * i + 1 ] = p2 - idOffset;
1506 
1507  }
1508 
1509  // done parsing edge section
1510  flag &= ~2;
1511  std::cout << "loaded." << std::endl;
1512  break;
1513  }
1514  }
1515 
1516  std::ifstream fstreamv ( fileName.c_str() );
1517  while ( getline ( fstreamv, line ) )
1518  {
1519 
1520  if ( line.find ( "volumeelements" ) != std::string::npos && flag & 8 )
1521  {
1522  getline ( fstreamv, line );
1523  std::stringstream parseLine ( line );
1524 
1525  parseLine >> numberVolumes;
1526  std::cout << "[readNetgenMesh] found " << numberVolumes << " volumes... " << std::flush;
1527 
1528  volumePointID.resize ( 4 * numberVolumes );
1529 
1530  for ( UInt i = 0; i < numberVolumes; i++ )
1531  {
1532  UInt matnr, np, p1, p2, p3, p4;
1533 
1534  getline ( fstreamv, line );
1535  std::stringstream parseLine ( line );
1536 
1537  parseLine >> matnr >> np >> p1 >> p2 >> p3 >> p4;
1538 
1539  volumePointID[ 4 * i ] = p1 - idOffset;
1540  volumePointID[ 4 * i + 1 ] = p2 - idOffset;
1541  volumePointID[ 4 * i + 2 ] = p3 - idOffset;
1542  volumePointID[ 4 * i + 3 ] = p4 - idOffset;
1543 
1544  ASSERT ( np == 4, "Error in readNetgenMesh: only tetrahedra elements supported" )
1545  }
1546  // done parsing volume section
1547  flag &= ~8;
1548  std::cout << "loaded." << std::endl;
1549  break;
1550  }
1551  }
1552  fstreamv.close();
1553 
1554 
1555  std::ifstream fstreamf ( fileName.c_str() );
1556  while ( getline ( fstreamf, line ) )
1557  {
1558 
1559  // TP 08/2008
1560  // surface elements section in a .vol file
1561  // is identified by key word "surfaceelements" in Netgen 4.5 (tested in RC2)
1562  // "surfaceelementsgi" is used in previous releases
1563  if ( ( line.find ("surfaceelements") != std::string::npos ) && flag & 4 )
1564  {
1565  getline ( fstreamf, line );
1566  std::stringstream parseLine ( line );
1567 
1568  parseLine >> numberBoundaryFaces;
1569 
1570  std::cout << "[readNetgenMesh] found " << numberBoundaryFaces
1571  << " boundary faces... " << std::flush;
1572 
1573  facePointID.resize ( 3 * numberBoundaryFaces );
1574  bcnsurf.resize ( numberBoundaryFaces + 1 );
1575  bcnsurf[ 0 ] = faceMarker.nullMarkerID();
1576 
1577 
1578  for ( UInt i = 0; i < numberBoundaryFaces; i++ )
1579  {
1580  UInt surfnr, bcnr, domin, domout, np, p1, p2, p3;
1581 
1582  getline ( fstreamf, line );
1583 
1584  if ( line.empty() )
1585  {
1586  // newer version of netgen inserts a blank line
1587  // after each line in "surface elements" section
1588  // make sure we ignore it
1589  // std::cout << "\nfound empty line in netgen file!" << std::endl;
1590  getline ( fstreamf, line );
1591  }
1592  std::stringstream parseLine ( line );
1593 
1594  parseLine >> surfnr >> bcnr >> domin >> domout >> np >> p1 >> p2 >> p3;
1595  p1 -= idOffset;
1596  p2 -= idOffset;
1597  p3 -= idOffset; //get the 0-based numbering
1598 
1599  facePointID[ 3 * i ] = p1;
1600  facePointID[ 3 * i + 1 ] = p2;
1601  facePointID[ 3 * i + 2 ] = p3;
1602 
1603  bcnsurf[ i + 1 ] = bcnr;
1604  // std::cout<<"[readNetgenMesh] bcnr = " << bcnr << std::endl;
1605  ASSERT ( np == 3, "Error in readNetgenMesh: only triangular surfaces supported" )
1606 
1607  //assume p1!=p2!=p3!=p1
1608  numberBoundaryVertices += boundaryPoint[ p1 ] ? 0 : 1;
1609  numberBoundaryVertices += boundaryPoint[ p2 ] ? 0 : 1;
1610  numberBoundaryVertices += boundaryPoint[ p3 ] ? 0 : 1;
1611  boundaryPoint[ p1 ] = boundaryPoint[ p2 ] = boundaryPoint[ p3 ] = true;
1612 
1613  /* here I set the boundary points marker
1614  note: this works only with my patch
1615  of strongerMarkerID
1616  Face flag is assigned to face points
1617  A point receives the "stronger" flag of the faces it belongs to
1618  */
1619  bcnpoints[ p1 ] = pointMarker.setStrongerMarkerID ( bcnpoints[ p1 ], bcnr );
1620  bcnpoints[ p2 ] = pointMarker.setStrongerMarkerID ( bcnpoints[ p2 ], bcnr );
1621  bcnpoints[ p3 ] = pointMarker.setStrongerMarkerID ( bcnpoints[ p3 ], bcnr );
1622 
1623  /* now I have the surface and points, so I can calculate
1624  the num of edges on the boundary, useful in case this is
1625  a quadratic tetra, I calculate the bcnr too using
1626  BareItemHandler<BareEdge> to store results
1627  */
1628  // the following is a complete list of edges in the 2D case
1629  // (only boundary edges in 3D)
1630  // may be useful even in the TWODIM case
1631  // (I've found this silly but easy way MM)
1632  BareEdge bed = setBareEdge ( p1, p2 );
1633 
1634  bihBedges.addIfNotThere ( bed, edgeMarker.nullMarkerID() );
1635  bihBedges[ bed ] = ( ID ) edgeMarker.setStrongerMarkerID ( bihBedges[ bed ], bcnr );
1636 
1637  bed = setBareEdge ( p2, p3 );
1638  bihBedges.addIfNotThere ( bed, edgeMarker.nullMarkerID() );
1639  bihBedges[ bed ] = ( ID ) edgeMarker.setStrongerMarkerID ( bihBedges[ bed ], bcnr );
1640 
1641  bed = setBareEdge ( p3, p1 );
1642  bihBedges.addIfNotThere ( bed, edgeMarker.nullMarkerID() );
1643  bihBedges[ bed ] = ( ID ) edgeMarker.setStrongerMarkerID ( bihBedges[ bed ], bcnr );
1644  }
1645  flag &= ~4;
1646  // in the 3D case the only way to know the number of edges on boundary faces
1647  // is to count them!
1648  numberBoundaryEdges = bihBedges.howMany();
1649  std::cout << "loaded." << std::endl;
1650  break;
1651  }
1652  }
1653 
1654  fstreamf.close();
1655 
1656  ASSERT ( flag == 0, "[readNetgenMesh] the mesh file does not have all the required sections." )
1657 
1658  std::cout << "[readNetgenMesh] computed " << numberBoundaryVertices << " boundary vertices" << std::endl;
1659 
1660  // Euler formulas
1661  numberFaces = 2 * numberVolumes + ( numberBoundaryFaces / 2 );
1662  numberEdges = numberVolumes + numberVertices + ( 3 * numberBoundaryFaces - 2 * numberBoundaryVertices ) / 4;
1663 
1664  // Be a little verbose
1665  if ( GeoShape::S_numPoints > 4 )
1666  {
1667  std::cout << "Quadratic Tetra Mesh (from Linear geometry)" << std::endl;
1668  numberPoints = numberVertices + numberEdges;
1669  numberBoundaryPoints = numberBoundaryVertices + numberBoundaryEdges; // I calculated the real numberBoundaryEdges before
1670  }
1671 
1672  else
1673  {
1674  std::cout << "Linear Tetra Mesh" << std::endl;
1675  numberPoints = numberVertices;
1676  numberBoundaryPoints = numberBoundaryVertices;
1677  }
1678 
1679  std::stringstream discardedLog;
1680  std::ostream& oStr = verbose ? std::cout : discardedLog;
1681  //points can be only vertices or on edges too
1682 
1683  std::cout << "Number of Vertices = " << std::setw ( 10 ) << numberVertices << std::endl
1684  << "Number of BVertices = " << std::setw ( 10 ) << numberBoundaryVertices << std::endl;
1685  oStr << "Number of Faces = " << std::setw ( 10 ) << numberFaces << std::endl
1686  << "Number of Boundary Faces = " << std::setw ( 10 ) << numberBoundaryFaces << std::endl
1687  << "Number of Edges = " << std::setw ( 10 ) << numberEdges << std::endl
1688  << "Number of Boundary Edges = " << std::setw ( 10 ) << numberBoundaryEdges << std::endl;
1689  std::cout << "Number of Points = " << std::setw ( 10 ) << numberPoints << std::endl
1690  << "Number of Boundary Points = " << std::setw ( 10 ) << numberBoundaryPoints << std::endl
1691  << "Number of Volumes = " << std::setw ( 10 ) << numberVolumes << std::endl;
1692 
1693  // Set all basic data structure
1694 
1695  // I store all Points
1696  mesh.setMaxNumPoints ( numberPoints, true );
1697  mesh.setMaxNumGlobalPoints ( numberPoints );
1698  mesh.setNumBPoints ( numberBoundaryPoints );
1699  mesh.setNumVertices ( numberVertices );
1700  mesh.setNumGlobalVertices ( numberVertices );
1701  mesh.setNumBVertices ( numberBoundaryVertices );
1702 
1703  // Only Boundary Edges
1704  mesh.setMaxNumEdges ( numberBoundaryEdges );
1705  mesh.setMaxNumGlobalEdges ( numberBoundaryEdges );
1706  mesh.setNumEdges ( numberEdges ); // Here the REAL number of edges (all of them)
1707  mesh.setNumBEdges ( numberBoundaryEdges );
1708 
1709  // Only Boundary Faces
1710  mesh.setMaxNumFaces ( numberBoundaryFaces );
1711  mesh.setMaxNumGlobalFaces ( numberBoundaryFaces );
1712  mesh.setNumFaces ( numberFaces ); // Here the REAL number of faces (all of them)
1713  mesh.setNumBFaces ( numberBoundaryFaces );
1714 
1715  mesh.setMaxNumVolumes ( numberVolumes, true );
1716  mesh.setMaxNumGlobalVolumes ( numberVolumes );
1717 
1718  mesh.setMarkerID ( regionFlag ); // Add Marker to list of Markers
1719 
1720  typename mesh_Type::point_Type* pointerPoint = 0;
1721  typename mesh_Type::edge_Type* pointerEdge = 0;
1722  typename mesh_Type::face_Type* pointerFace = 0;
1723  typename mesh_Type::volume_Type* pointerVolume = 0;
1724 
1725  // addPoint()/Face()/Edge() returns a reference to the last stored point
1726  // I use that information to set all point info, by using a pointer.
1727 
1728  std::cout << "[importerMesh3D] boundaryPoint.size() = " << boundaryPoint.size()
1729  << ", bcnpoints.size() = " << bcnpoints.size()
1730  << std::endl;
1731 
1732  for ( UInt i = 0; i < numberVertices; i++ )
1733  {
1734  pointerPoint = &mesh.addPoint ( boundaryPoint[ i ], true ); //true if boundary point
1735 
1736  pointerPoint->setId ( i );
1737  pointerPoint->setMarkerID ( bcnpoints[ i ] );
1738  pointerPoint->x() = pointCoordinates[ nDimensions * i ];
1739  pointerPoint->y() = pointCoordinates[ nDimensions * i + 1 ];
1740  pointerPoint->z() = pointCoordinates[ nDimensions * i + 2 ];
1741  }
1742 
1743  std::cout << "[importerMesh3D] added points." << std::endl;
1744 
1745  /*
1746  here I set the real boundary edges that I stored
1747  in bihBedges
1748  */
1749 
1750  MeshElementBareHandler<BareEdge>::const_iterator bedge = bihBedges.begin();
1751 
1752  for ( UInt i = 0; i < numberBoundaryEdges; i++ )
1753  {
1754  UInt p1, p2;
1755 
1756  pointerEdge = &mesh.addEdge ( true ); // Only boundary edges.
1757  pointerEdge->setId ( i );
1758  pointerEdge->setMarkerID ( markerID_Type ( bedge->second ) );
1759  p1 = bedge->first.first;
1760  p2 = bedge->first.second;
1761  pointerEdge->setPoint ( 0, mesh.point ( p1 ) ); // set edge conn.
1762  pointerEdge->setPoint ( 1, mesh.point ( p2 ) ); // set edge conn.
1763 
1764  bedge++;
1765  }
1766 
1767  std::cout << "[importerMesh3D] added edges." << std::endl;
1768 
1769  for ( UInt i = 0; i < numberVolumes; i++ )
1770  {
1771  UInt p1, p2, p3, p4;
1772 
1773  pointerVolume = &mesh.addVolume();
1774  pointerVolume->setId ( i );
1775  p1 = volumePointID[ 4 * i ];
1776  p2 = volumePointID[ 4 * i + 1 ];
1777  p3 = volumePointID[ 4 * i + 2 ];
1778  p4 = volumePointID[ 4 * i + 3 ];
1779  pointerVolume->setPoint ( 0, mesh.point ( p1 ) );
1780  pointerVolume->setPoint ( 1, mesh.point ( p2 ) );
1781  pointerVolume->setPoint ( 2, mesh.point ( p3 ) );
1782  pointerVolume->setPoint ( 3, mesh.point ( p4 ) );
1783  }
1784 
1785  std::cout << "[importerMesh3D] added volumes." << std::endl;
1786 
1787  for ( UInt i = 0; i < numberBoundaryFaces; i++ )
1788  {
1789  UInt p1, p2, p3;
1790 
1791  pointerFace = &mesh.addFace ( true ); // Only boundary faces
1792  p1 = facePointID[ 3 * i ];
1793  p2 = facePointID[ 3 * i + 1 ];
1794  p3 = facePointID[ 3 * i + 2 ];
1795 
1796  pointerFace->setMarkerID ( markerID_Type ( bcnsurf[ i + 1 ] ) );
1797  pointerEdge->setId ( i );
1798  pointerFace->setPoint ( 0, mesh.point ( p1 ) ); // set face conn.
1799  pointerFace->setPoint ( 1, mesh.point ( p2 ) ); // set face conn.
1800  pointerFace->setPoint ( 2, mesh.point ( p3 ) ); // set face conn.
1801  }
1802  std::cout << "[importerMesh3D] added faces." << std::endl;
1803 
1804  // This part is to build a P2 mesh from a P1 geometry
1805 
1806  if ( GeoShape::S_numPoints > 4 )
1807  {
1808  MeshUtility::p2MeshFromP1Data ( mesh );
1809  }
1810 
1811  // Test mesh
1812  Switch sw;
1813 
1814  ///// CORRECTION JFG
1815  //if (mesh.check(1, true,true))done=0;
1816 
1817  if ( !checkMesh3D ( mesh, sw, true, verbose, oStr, std::cerr, oStr ) )
1818  {
1819  std::abort(); // CORRECTION JFG
1820  }
1821 
1822  Real vols[ 3 ];
1823 
1824  getVolumeFromFaces ( mesh, vols, oStr );
1825 
1826  oStr << "Volume enclosed by the mesh computed by integration on boundary faces" << std::endl
1827  << "INT(X) INT(Y) INT(Z) <- they should be equal and equal to" << std::endl
1828  << " the volume enclosed by the mesh " << std::endl
1829  << vols[ 0 ] << " " << vols[ 1 ] << " " << vols[ 2 ] << std::endl
1830  << "Boundary faces are defining aclosed surface if "
1831  << testClosedDomain ( mesh, oStr ) << std::endl
1832  << " is (almost) zero" << std::endl;
1833 
1834  return true;
1835 } // Function readNetgenMesh
1836 
1837 //! saveNetgenSolution -
1838 /*!
1839  Ripped "from src/ng431/libsrc/interface/importsolution.cpp"
1840 
1841  @param fileName, the name of the mesh file to read.
1842  @param U,
1843  @param fctname, default to "u".
1844  @return void.
1845 */
1846 
1847 template <typename VectorType>
1848 void
1849 saveNetgenSolution (std::string fileName,
1850  const VectorType& solution,
1851  std::string functionName = "u")
1852 {
1853  std::ofstream of (fileName.c_str() );
1854 
1855  of << "solution " << functionName << " -size =" << solution.size() << " -components=1 -type=nodal" << std::endl;
1856 
1857  for ( UInt i = 0; i < solution.size(); ++i )
1858  {
1859  of << solution ( i ) << std::endl;
1860  }
1861 
1862  of.close();
1863 } // Function saveNetgenSolution
1864 
1865 } // Namespace LifeV
1866 
1867 #endif /* IMPORTERMESH3D_H */
ASSERT_PRE0
#define ASSERT_PRE0(X, A)
Definition: LifeAssert.hpp:72
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers
Definition: ImporterMesh3D.hpp:72
LifeV::Switch
I use standard constructor/destructors.
Definition: Switch.hpp:50
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers::i1
UInt i1
Definition: ImporterMesh3D.hpp:75
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FaceHelp::ibc
Int ibc
Definition: ImporterMesh3D.hpp:83
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FaceHelp
Definition: ImporterMesh3D.hpp:79
LifeV::Int
int32_type Int
Generic integer data.
Definition: LifeV.hpp:188
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers::i4
UInt i4
Definition: ImporterMesh3D.hpp:75
LifeV::markerID_Type
ID markerID_Type
markerID_Type is the type used to store the geometric entity marker IDs
Definition: Marker.hpp:81
ERROR_MSG
#define ERROR_MSG(A)
Definition: LifeAssert.hpp:69
LifeV::BareEdge
The Edge basis class.
Definition: MeshElementBare.hpp:176
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers::ibc
Int ibc
Definition: ImporterMesh3D.hpp:76
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FaceHelp::i
UInt i[9]
Definition: ImporterMesh3D.hpp:82
LifeV::readMppFileHead
bool readMppFileHead(std::ifstream &myStream, UInt &numberVertices, UInt &numberBoundaryVertices, UInt &numberBoundaryFaces, UInt &numberBoundaryEdges, UInt &numberVolumes)
readMppFileHead - reads mesh++ Tetra meshes.
Definition: ImporterMesh3D.cpp:50
LifeV::MeshElementBareHandler::containerConstIterator_Type
container_Type::const_iterator containerConstIterator_Type
Definition: MeshElementBare.hpp:565
LifeV::Assert::error
Assert & error(const char *strMsg=0)
Definition: LifeAssertSmart.hpp:326
LifeV::Vector
boost::numeric::ublas::vector< Real > Vector
Definition: ImporterMesh3D.hpp:66
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers::i3
UInt i3
Definition: ImporterMesh3D.hpp:75
LifeV::anonymous_namespace{ImporterMesh3D.hpp}::FiveNumbers::i2
UInt i2
Definition: ImporterMesh3D.hpp:75
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191