readMesh3D.hpp 52.4 KB
Newer Older
prudhomm's avatar
prudhomm committed
1
/*
2
  This file is part of the LifeV library
3
  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politecnico di Milano
4
5
6
7
8
9
10
11
12
13
14
15
16
17

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
prudhomm's avatar
prudhomm committed
18
*/
prudhomm's avatar
prudhomm committed
19
20
#ifndef _READMESH3D_HH_
#define _READMESH3D_HH_
21
22
23
24
25
26
27
28

#include <sstream>
#include <algorithm>

#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/lambda/if.hpp>

prudhomm's avatar
prudhomm committed
29
#include <life/lifecore/debug.hpp>
30
31


prudhomm's avatar
prudhomm committed
32
33
34
#include <life/lifemesh/regionMesh3D.hpp>
#include <life/lifecore/util_string.hpp>
#include <life/lifefilters/mesh_util.hpp>
35
#include <life/lifefilters/selectMarker.hpp>
prudhomm's avatar
prudhomm committed
36
37
38

namespace LifeV
{
39
/*----------------------------------------------------------------------*
40
41
42
43
44
45
46
  |
  | Added support for numbering from 1
  ! Added markers
  |
  | #Mesh readers
  |
  *----------------------------------------------------------------------*/
47
48
49
50
/********************************************************************************
MESH BUILDERS
********************************************************************************/

51
52
53
54
55
// template<typename T>
// T Max(T & i1, T & i2,T & i3)
// {
//     return std::max(std::max(i1,i2),i3);
// }
56

57
58
59
60
61
// template<typename T>
// T Max(T & i1, T & i2, T & i3, T & i4)
// {
//     return std::max(Max(i1,i2,i3),i4);
// }
62

63
bool
64
65
66
67
68
69
readMppFileHead( std::ifstream & mystream,
                 UInt & numVertices,
                 UInt & numBVertices,
                 UInt & numBFaces,
                 UInt & numBEdges,
                 UInt & numVolumes );
70
71
72
73
74
75

//
//=================================================================
// ReadmppFile It reads mesh++ Tetra meshes. It convert them into
// Quadratic Tetra if needed it return false if  mesh check is uncessessfull
//
76
template <typename GeoShape, typename MC>
prudhomm's avatar
prudhomm committed
77
bool
78
79
80
readMppFile( RegionMesh3D<GeoShape, MC> & mesh,
             const std::string & filename,
             EntityFlag regionFlag,
81
             bool verbose = false )
82
83
84
85
86
87
88
{
    unsigned done = 0;
    std::string line;
    Real x, y, z;
    int ity, ity_id;
    UInt p1, p2, p3, p4;
    UInt nVe( 0 ), nBVe( 0 ), nFa( 0 ), nBFa( 0 ), nPo( 0 ), nBPo( 0 );
simone's avatar
simone committed
89
90
    UInt nVo( 0 ), nBEd( 0 );
    UInt nEd( 0 );
91
    UInt i;
winkelma's avatar
winkelma committed
92

93
94
95
    std::stringstream discardedLog;
    std::ostream& oStr = verbose ? std::cout : discardedLog;

96
97
98
99
100
101
102
103
104
105
106
    ASSERT_PRE0( GeoShape::Shape == TETRA , "ReadMppFiles reads only tetra meshes" ) ;

    ASSERT_PRE0( GeoShape::Shape == TETRA, "Sorry, ReadMppFiles reads only tetra meshes" );

    ASSERT_PRE0( GeoShape::numVertices <= 6, "Sorry, ReadMppFiles handles only liner&quad tetras" );

    // open stream to read header

    std::ifstream hstream( filename.c_str() );
    if ( hstream.fail() )
    {
107
108
        std::cerr << " Error in readMpp: File " << filename
                  << " not found or locked" << std::endl;
109
110
        abort();
    }
111
    std::cout << "Reading Mesh++ file" << std::endl;
112
113
114
115
116
117
118
119
120
121
122
    if ( ! readMppFileHead( hstream, nVe, nBVe, nBFa, nBEd, nVo ) )
    {
        std::cerr << " Error While reading mesh++ file headers" << std::endl;
        ABORT() ;
    }
    hstream.close();

    //Reopen the stream: I know it is stupid but this is how it goes
    std::ifstream mystream( filename.c_str() );
    if ( mystream.fail() )
    {
123
124
        std::cerr << " Error in readMpp: File " << filename
                  << " not found or locked" << std::endl;
125
126
127
128
129
130
131
132
133
134
135
        abort();
    }

    // Euler formulas
    nFa = 2 * nVo + ( nBFa / 2 );
    nEd = nVo + nVe + ( 3 * nBFa - 2 * nBVe ) / 4;

    // Be a little verbose
    if ( GeoShape::numPoints > 4 )
    {

136
137
        std::cout << "Quadratic Tetra  Mesh (from Linear geometry)"
                  << std::endl;
138
139
140
141
142
143
144
145
146
        nPo = nVe + nEd;
        nBPo = nBVe + nBEd;
    }
    else
    {
        std::cout << "Linear Tetra Mesh" << std::endl;
        nPo = nVe;
        nBPo = nBVe;
    }
147
148
149
150
151
152
153
154
155
    std::cout << "#Vertices = "          << std::setw(10) << nVe
              << "  #BVertices       = " << std::setw(10) << nBVe << std::endl;
    oStr      << "#Faces    = "          << std::setw(10) << nFa
              << "  #Boundary Faces  = " << std::setw(10) << nBFa << std::endl;
    oStr      << "#Edges    = "          << std::setw(10) << nEd
              << "  #Boundary Edges  = " << std::setw(10) << nBEd << std::endl;
    std::cout << "#Points   = "          << std::setw(10) << nPo
              << "  #Boundary Points = " << std::setw(10) << nBPo << std::endl;
    std::cout << "#Volumes  = "          << std::setw(10) << nVo  << std::endl;
156
157
158
159
160

    // Set all basic data structure

    // I store all Points
    mesh.setMaxNumPoints( nPo, true );
simone's avatar
simone committed
161
162
163
164
    mesh.setMaxNumGlobalPoints( nPo );
    mesh.setNumBPoints  ( nBPo );
    mesh.setNumVertices ( nVe );
    mesh.setNumGlobalVertices(nVe);
simone's avatar
simone committed
165
    mesh.setNumBVertices( nBVe );
166
    // Only Boundary Edges (in a next version I will allow for different choices)
simone's avatar
simone committed
167
168
169
170
    mesh.setMaxNumEdges ( nEd );
    mesh.setMaxNumGlobalEdges ( nEd );
    mesh.setNumEdges    ( nEd ); // Here the REAL number of edges (all of them)
    mesh.setNumBEdges   ( nBEd );
171
    // Only Boundary Faces
simone's avatar
simone committed
172
173
174
175
    mesh.setMaxNumFaces ( nBFa );
    mesh.setMaxNumGlobalFaces ( nBFa );
    mesh.setNumFaces    ( nFa ); // Here the REAL number of edges (all of them)
    mesh.setNumBFaces   ( nBFa );
176
177

    mesh.setMaxNumVolumes( nVo, true );
simone's avatar
simone committed
178
    mesh.setMaxNumGlobalVolumes( nVo);
179
180
181

    mesh.setMarker( regionFlag ); // Mark the region

simone's avatar
simone committed
182

183
184
185
186
187
188
189
190
191
192
193
194
    typename RegionMesh3D<GeoShape, MC>::PointType * pp = 0;
    typename RegionMesh3D<GeoShape, MC>::EdgeType * pe = 0;
    typename RegionMesh3D<GeoShape, MC>::FaceType * pf = 0;
    typename RegionMesh3D<GeoShape, MC>::VolumeType * pv = 0;
    // addPoint()/Face()/Edge() returns a reference to the last stored point
    // I use that information to set all point info, by using a pointer.
    UInt count = 0;
    long int ibc;
    while ( next_good_line( mystream, line ).good() )
    {
        if ( line.find( "odes" ) != std::string::npos )
        {
195

196
197
198
199
            std::string node_s = line.substr( line.find_last_of( ":" ) + 1 );
            //      _numVertices=atoi(node_s);
            for ( i = 0;i < nVe;i++ )
            {
200
#ifdef OLDMPPFILE
201
                mystream >> x >> y >> z >> ity >> ibc;
202
#else
203
204
205
206

                mystream >> x >> y >> z >> ity >> ity_id;
                if ( ity != 3 )
                    mystream >> ibc;
207
#endif
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222

                if ( ity != 3 )
                {
                    ++count;
                    pp = &mesh.addPoint( true ); // Boundary point. Boundary switch set by the mesh method.
                    pp->setMarker( EntityFlag( ibc ) );
                }
                else
                {
                    pp = &mesh.addPoint( false );
                }
                pp->x() = x;
                pp->y() = y;
                pp->z() = z;
                pp->setMarker( EntityFlag( ibc ) );
223
224
225
226
227
228

                pp->setId     ( i + 1 );
                pp->setLocalId( i + 1 );

                mesh.localToGlobalNode().insert(std::make_pair(i+1, i+1));
                mesh.globalToLocalNode().insert(std::make_pair(i+1, i+1));
229
            }
230
            oStr << "Vertices Read " << std::endl;
231
232
233
234
235
236
            done++;
            if ( count != nBVe )
                std::cerr << "NumB points inconsistent !" << std::endl;
        }
        if ( line.find( "iangular" ) != std::string::npos )
        {
237
            oStr << "Reading Bfaces " << std::endl;
238
239
240
241
            std::string node_s = line.substr( line.find_last_of( ":" ) + 1 );
            // _numBFaces=atoi(node_s);
            for ( i = 0;i < nBFa;i++ )
            {
242
#ifdef OLDMPPFILE
243
                mystream >> p1 >> p2 >> p3 >> ity >> ibc;
244
#else
245
246

                mystream >> p1 >> p2 >> p3 >> ity >> ity_id >> ibc;
247
#endif
248
249
250
251
252
253
254
255

                pf = &( mesh.addFace( true ) ); // Only boundary faces

                pf->setMarker( EntityFlag( ibc ) );
                pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
            }
256
            oStr << "Boundary Faces Read " << std::endl;
257
258
259
260
            done++;
        }
        if ( line.find( "Sides" ) != std::string::npos )
        {
261
            oStr << "Reading Bedges " << std::endl;
262
263
264
265
            std::string node_s = line.substr( line.find_last_of( ":" ) + 1 );
            //_numBEdges=atoi(node_s);
            for ( i = 0;i < nBEd;i++ )
            {
266
#ifdef OLDMPPFILE
267
                mystream >> p1 >> p2 >> ity >> ibc;
268
#else
269
270

                mystream >> p1 >> p2 >> ity >> ity_id >> ibc;
271
#endif
prudhomm's avatar
prudhomm committed
272

273
274
275
276
277
                pe = &mesh.addEdge( true ); // Only boundary edges.
                pe->setMarker( EntityFlag( ibc ) );
                pe->setPoint( 1, mesh.point( p1 ) ); // set edge conn.
                pe->setPoint( 2, mesh.point( p2 ) ); // set edge conn.
            }
278
            oStr << "Boundary Edges Read " << std::endl;
279
280
281
282
283
            done++;
        }
        count = 0;
        if ( line.find( "etrahedral" ) != std::string::npos )
        {
284
            oStr << "Reading Volumes " << std::endl;
285
286
287
288
289
            std::string node_s = line.substr( line.find_last_of( ":" ) + 1 );
            for ( i = 0; i < nVo; i++ )
            {
                mystream >> p1 >> p2 >> p3 >> p4;
                pv = &mesh.addVolume();
simone's avatar
simone committed
290
291
292
                pv->setId     ( i + 1 );
                pv->setLocalId( i + 1);
//                pv->id() = i + 1;
293
294
295
296
297
298
                pv->setPoint( 1, mesh.point( p1 ) );
                pv->setPoint( 2, mesh.point( p2 ) );
                pv->setPoint( 3, mesh.point( p3 ) );
                pv->setPoint( 4, mesh.point( p4 ) );
                count++;
            }
299
            oStr << count << " Volume elements Read" << std::endl;
300
301
302
303
304
            done++;
        }

    }
    // This part is to build a P2 mesh from a P1 geometry
prudhomm's avatar
prudhomm committed
305

306
307
308
    if ( GeoShape::numPoints > 4 )
        p1top2( mesh );
    mystream.close();
prudhomm's avatar
prudhomm committed
309

310
311
    // Test mesh
    Switch sw;
prudhomm's avatar
prudhomm committed
312

313
314
    ///// CORRECTION JFG
    //if (mesh.check(1, true,true))done=0;
prudhomm's avatar
prudhomm committed
315

316
    if ( !checkMesh3D( mesh, sw, true, verbose, oStr, std::cerr, oStr ) )
317
        abort(); // CORRECTION JFG
318

319
    Real vols[ 3 ];
320
321
    getVolumeFromFaces( mesh, vols, oStr );
    oStr << "   VOLUME ENCLOSED BY THE MESH COMPUTED BY INTEGRATION ON" <<
322
        " BOUNDARY FACES" << std::endl;
323
    oStr << "INT(X)     INT(Y)      INT(Z) <- they should be equal and equal to" << std::endl <<
324
        "                                 the voulume enclosed by the mesh " << std::endl;
325
    oStr << vols[ 0 ] << " " << vols[ 1 ] << " " << vols[ 2 ] << std::endl;
prudhomm's avatar
prudhomm committed
326

327
    oStr << "   BOUNDARY FACES ARE DEFINING A CLOSED SURFACE IF " << testClosedDomain( mesh, oStr ) << std::endl <<
328
        " IS (ALMOST) ZERO" << std::endl;
prudhomm's avatar
prudhomm committed
329

330
    return done == 4 ;
prudhomm's avatar
prudhomm committed
331

332
333
334
}

/*
335
336
337
  INRIA MESH FILE READERS
  29/06/2002 L.F.
*/
338
339
//! INRIAMesh used either spaces or CR as separators
/*! It sucks, but this is the way it is! This little function should help handling it. It gets an
340
  integer field from the std::string line if it is not empty, otherwise from the input stream.
341

342
  It assumes that the std::string is either empty or it contains and integer!!! No check is made to verify this.
343
344
*/

345
int nextIntINRIAMeshField( std::string const & line, std::istream & mystream );
346
347

bool
348
349
350
351
352
353
readINRIAMeshFileHead( std::ifstream & mystream,
                       UInt & numVertices,
                       UInt & numBVertices,
                       UInt & numBFaces,
                       UInt & numBEdges,
                       UInt & numVolumes,
354
                       UInt & numStoredFaces,
355
                       ReferenceShapes & shape,
356
357
                       InternalEntitySelector
                       iSelect=InternalEntitySelector());
358

359
360
/*!
  read an INRIA mesh
361
*/
362
363
364
struct FiveNumbers
{
public:
365
366
    UInt i1,i2,i3,i4;
    long int ibc;
367
};
368
template <typename GeoShape, typename MC>
prudhomm's avatar
prudhomm committed
369
bool
370
371
372
373
374
readINRIAMeshFile( RegionMesh3D<GeoShape, MC>&      mesh,
                   std::string const&               filename,
                   EntityFlag                       regionFlag,
                   bool verbose                   = false,
                   InternalEntitySelector iSelect = InternalEntitySelector())
375
376
377
378
379
{
    unsigned done = 0;
    std::string line;
    Real x, y, z;
    UInt p1, p2, p3, p4, p5, p6, p7, p8;
simone's avatar
simone committed
380
381
    UInt nVe( 0 ), nFa( 0 ), nBFa( 0 ), nPo( 0 ), nBPo( 0 );
    UInt nBVe(0);
382
    UInt numStoredFaces(0);
simone's avatar
simone committed
383
384
    UInt nVo( 0 ), nBEd( 0 );
    UInt nEd;
385
    UInt i;
386
    ReferenceShapes shape(NONE);
387
388
    std::vector<FiveNumbers> faceHelp;
    typename std::vector<FiveNumbers>::iterator faceHelpIterator;
389
390
391
    std::stringstream discardedLog;
    std::ostream& oStr = verbose ? std::cout : discardedLog;

392
393
394
    // open stream to read header

    std::ifstream hstream( filename.c_str() );
395
396
397
398
    if (verbose)
        {
            std::cout<<"Reading form file "<<filename<< std::endl;
        }
399

400
401
    if ( hstream.fail() )
    {
402
403
404
        std::cerr << " Error in readINRIAMeshFile: File " << filename
                  << " not found or locked"
                  << std::endl;
405
        abort();
406
    }
407
    if (verbose) std::cout << "Reading INRIA mesh file" << filename << std::endl;
lformaggia's avatar
lformaggia committed
408
    if ( ! readINRIAMeshFileHead( hstream, nVe, nBVe, nBFa, nBEd, nVo, numStoredFaces,shape,
409
                                  iSelect) )
410
411
412
    {
        std::cerr << " Error While reading INRIA mesh file headers" << std::endl;
        ABORT() ;
413
    }
414
    hstream.close();
lformaggia's avatar
lformaggia committed
415

416
417
418
419
    //Reopen the stream: I know it is stupid but this is how it goes
    std::ifstream mystream( filename.c_str() );
    if ( mystream.fail() )
    {
420
421
        std::cerr << " Error in readINRIAMeshFile: File " << filename
                  << " not found or locked" << std::endl;
422
        abort();
423
424
    }

425
426
427
428
    ASSERT_PRE0( GeoShape::Shape == shape, "INRIA Mesh file and mesh element shape is not consistent" );

    // Euler formulas to get number of faces and number of edges
    nFa = 2 * nVo + ( nBFa / 2 );
simone's avatar
simone committed
429
430
431
432
433
434
435
    int num1  = nVe + nVo;
    int num2  = nBVe;
    int num3  = nBFa;

    nEd = (3*num3 - 2*num2)/4 + num1;

//    nEd = (int) nVo + nVe + ( 3 * nBFa + dummy ) / 4;
436
437
438
439
440

    // Be a little verbose
    switch ( shape )
    {

441
442
443
        case HEXA:
            ASSERT_PRE0( GeoShape::numPoints == 8, "Sorry I can read only bilinear Hexa meshes" );
            std::cout << "Linear Hexa Mesh" << std::endl;
444
445
            nPo = nVe;
            nBPo = nBVe;
446
447
448
449
450
451
452
453
454
            break;
        case TETRA:
            if ( GeoShape::numPoints > 4 )
            {
                //    if (GeoShape::numPoints ==6 ){
                std::cout << "Quadratic Tetra  Mesh (from Linear geometry)" << std::endl;
                nPo = nVe + nEd;
                // nBPo=nBVe+nBEd; // FALSE : nBEd is not known at this stage in a INRIA file (JFG 07/2002)
                // I use the relation  nBVe + nBFa - 2 = nBEd, But, is it general (hole...) ???? (JFG 07/2002)
455
456
                nBEd=(int( nBVe + nBFa - int(2) )>0?( nBVe + nBFa - 2 ):0);
                nBPo = (int(nBVe + ( nBVe + nBFa - int(2) ))>0?nBVe + ( nBVe + nBFa - 2 ):0);
457
458
459
            }
            else
            {
460
461
462
                if (verbose)
                    std::cout << "Linear Tetra Mesh" << std::endl;

463
464
                nPo = nVe;
                nBPo = nBVe;
465
                nBEd=(int( nBVe + nBFa - int(2) )>0?( nBVe + nBFa - 2 ):0);
466
467
468
469
            }
            break;
        default:
            ERROR_MSG( "Current version of INRIA Mesh file reader only accepts TETRA and HEXA" );
470
    }
prudhomm's avatar
prudhomm committed
471

472
    oStr << "#Vertices = "          << std::setw(10) << nVe
473
              << "  #BVertices       = " << std::setw(10) << nBVe << std::endl;
474
    oStr << "#Faces    = "          << std::setw(10) << nFa
475
476
              << "  #Boundary Faces  = " << std::setw(10) << nBFa << std::endl
              << "#Stored Faces = " << std::setw(10) << numStoredFaces<< std::endl;
477
    oStr << "#Edges    = "          << std::setw(10) << nEd
478
              << "  #Boundary Edges  = " << std::setw(10) << nBEd << std::endl;
479
    oStr << "#Points   = "          << std::setw(10) << nPo
480
              << "  #Boundary Points = " << std::setw(10) << nBPo << std::endl;
481
    oStr << "#Volumes  = "          << std::setw(10) << nVo  << std::endl;
482
483
484
485

    // Set all basic data structure

    // I store all Points
simone's avatar
simone committed
486
487
488
489
490
491
    mesh.setMaxNumPoints   ( nPo, true );
    mesh.setMaxNumGlobalPoints( nPo );
    mesh.setNumBPoints     ( nBPo );
    mesh.setNumVertices    ( nVe );
    mesh.setNumGlobalVertices(nVe);
    mesh.setNumBVertices   ( nBVe );
492
    // Only Boundary Edges (in a next version I will allow for different choices)
simone's avatar
simone committed
493
494
495
496
    mesh.setMaxNumEdges    ( nBEd );
    mesh.setMaxNumGlobalEdges ( nEd );
    mesh.setNumEdges       ( nEd ); // Here the REAL number of edges (all of them)
    mesh.setNumBEdges      ( nBEd );
497
    // Only Boundary Faces
simone's avatar
simone committed
498
499
500
501
    mesh.setMaxNumFaces    ( numStoredFaces );
    mesh.setMaxNumGlobalFaces ( nBFa );
    mesh.setNumFaces       ( nFa ); // Here the REAL number of faces (all of them)
    mesh.setNumBFaces      ( nBFa );
502

simone's avatar
simone committed
503
504
    mesh.setMaxNumVolumes  ( nVo, true );
    mesh.setMaxNumGlobalVolumes( nVo);
505

simone's avatar
simone committed
506
    mesh.setMarker         ( regionFlag ); // Add Marker to list of Markers
507

508
509
510
511
    typedef typename RegionMesh3D<GeoShape, MC>::PointType  PointType;
    typedef typename RegionMesh3D<GeoShape, MC>::VolumeType VolumeType;


512
513
514
515
516
517
518
519
    typename RegionMesh3D<GeoShape, MC>::PointType * pp = 0;
    typename RegionMesh3D<GeoShape, MC>::EdgeType * pe = 0;
    typename RegionMesh3D<GeoShape, MC>::FaceType * pf = 0;
    typename RegionMesh3D<GeoShape, MC>::VolumeType * pv = 0;
    // addPoint()/Face()/Edge() returns a reference to the last stored point
    // I use that information to set all point info, by using a pointer.
    UInt count = 0;
    long int ibc;
520

521
522
    // To account for internal faces
    if (numStoredFaces > nBFa){
523
524
525
526
        faceHelp.resize(numStoredFaces - nBFa);
        faceHelpIterator=faceHelp.begin();
        oStr<<"WARNING: The mesh file (apparently) contains "<<numStoredFaces - nBFa<<" internal faces"<<std::endl;

527
    }
528

529
530
531
532
533
    while ( next_good_line( mystream, line ).good() )
    {
        if ( line.find( "Vertices" ) != std::string::npos )
        {
            nextIntINRIAMeshField( line.substr( line.find_last_of( "s" ) + 1 ), mystream );
simone's avatar
simone committed
534
            for ( i = 0; i < nVe; i++ )
535
            {
536
                mystream >> x >> y >> z >> ibc;
simone's avatar
simone committed
537
538
539

//                if (ibc == 1 ) ibc = 100;

lformaggia's avatar
lformaggia committed
540
                if ( !iSelect(EntityFlag(ibc)))
541
                {
542
543
544
                    ++count;
                    pp = &mesh.addPoint( true ); // Boundary point. Boundary switch set by the mesh method.
                    pp->setMarker( EntityFlag( ibc ) );
545
                }
546
                else
547
                {
lformaggia's avatar
lformaggia committed
548
                    pp = &mesh.addPoint( false );
549
                }
simone's avatar
simone committed
550
551
                pp->setId     ( i + 1 );
                pp->setLocalId( i + 1 );
552
553
554
555
                pp->x() = x;
                pp->y() = y;
                pp->z() = z;
                pp->setMarker( EntityFlag( ibc ) );
simone's avatar
simone committed
556

557
558
                mesh.localToGlobalNode().insert(std::make_pair(i+1, i+1));
                mesh.globalToLocalNode().insert(std::make_pair(i+1, i+1));
559
            }
560
            oStr << "Vertices Read " << std::endl;
561
            oStr << "size of the node storage is " << count*sizeof(PointType)/1024./1024. << std::endl;
562
563
            done++;
            if ( count != nBVe )
lformaggia's avatar
lformaggia committed
564
                std::cerr << "NumB points inconsistent !" << std::endl;
565
        }
lformaggia's avatar
lformaggia committed
566

567
        if ( line.find( "Triangles" ) != std::string::npos ){
568
569
570
            nextIntINRIAMeshField( line.substr( line.find_last_of( "s" ) + 1 ), mystream );
            oStr << "Reading Bfaces " << std::endl;
            for ( i = 0;i < numStoredFaces;i++ )
lformaggia's avatar
lformaggia committed
571
            {
572
                mystream >> p1 >> p2 >> p3 >> ibc;
lformaggia's avatar
lformaggia committed
573

574
                if (numStoredFaces > nBFa){
lformaggia's avatar
lformaggia committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
                    if (mesh.point( p1 ).boundary()&&mesh.point( p2 ).boundary()&&
                        mesh.point( p3 ).boundary()){
                        pf = &( mesh.addFace( true ) ); // Boundary faces
                        pf->setMarker( EntityFlag( ibc ) );
                        pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                        pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                        pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.

                    } else {
                        faceHelpIterator->i1=p1;
                        faceHelpIterator->i2=p2;
                        faceHelpIterator->i3=p3;
                        faceHelpIterator->ibc=ibc;
                        ++faceHelpIterator;
                    }
590
                } else {
lformaggia's avatar
lformaggia committed
591
592
593
594
595
596

                    pf = &( mesh.addFace( true ) ); // Only boundary faces
                    pf->setMarker( EntityFlag( ibc ) );
                    pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                    pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                    pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
597
598
599
600
601
602
603
604
605
                }
            }
            for (faceHelpIterator=faceHelp.begin();faceHelpIterator!=faceHelp.end();
                 ++faceHelpIterator){
                p1=faceHelpIterator->i1;
                p2=faceHelpIterator->i2;
                p3=faceHelpIterator->i3;
                ibc=faceHelpIterator->ibc;
                pf = &( mesh.addFace( false ) ); // INTERNAL FACE
606
607
608
609
                pf->setMarker( EntityFlag( ibc ) );
                pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
610
611
612
613
614
615
            }

            oStr << "Boundary Faces Read " << std::endl;
            done++;
        }

lformaggia's avatar
lformaggia committed
616
        if ( line.find( "Quadrilaterals" ) != std::string::npos ){
617
618
619
            nextIntINRIAMeshField( line.substr( line.find_last_of( "s" ) + 1 ), mystream );
            oStr << "Reading Bfaces " << std::endl;
            for ( i = 0;i < nBFa;i++ )
620
            {
621
622
623
                mystream >> p1 >> p2 >> p3 >> p4 >> ibc;

                if (numStoredFaces > nBFa){
lformaggia's avatar
lformaggia committed
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
                    if (mesh.point( p1 ).boundary()&&mesh.point( p2 ).boundary()&&
                        mesh.point( p3 ).boundary()){
                        pf = &( mesh.addFace( true ) ); // Boundary faces
                        pf->setMarker( EntityFlag( ibc ) );
                        pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                        pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                        pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
                        pf->setPoint( 4, mesh.point( p4 ) ); // set face conn.

                    } else {
                        faceHelpIterator->i1=p1;
                        faceHelpIterator->i2=p2;
                        faceHelpIterator->i3=p3;
                        faceHelpIterator->i4=p4;
                        faceHelpIterator->ibc=ibc;
                        ++faceHelpIterator;
                    }
641
                } else {
lformaggia's avatar
lformaggia committed
642
643
644
645
646
647
                    pf = &( mesh.addFace( true ) ); // Only boundary faces
                    pf->setMarker( EntityFlag( ibc ) );
                    pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                    pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                    pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
                    pf->setPoint( 4, mesh.point( p4 ) ); // set face conn.
648
649
650
651
652
653
654
655
656
657
658
                }
            }
            oStr << "Boundary Faces Read " << std::endl;
            for (faceHelpIterator=faceHelp.begin();faceHelpIterator!=faceHelp.end();
                 ++faceHelpIterator){
                p1=faceHelpIterator->i1;
                p2=faceHelpIterator->i2;
                p3=faceHelpIterator->i3;
                p4=faceHelpIterator->i4;
                ibc=faceHelpIterator->ibc;
                pf = &( mesh.addFace( false ) ); // INTERNAL FACE
659
660
661
662
663
                pf->setMarker( EntityFlag( ibc ) );
                pf->setPoint( 1, mesh.point( p1 ) ); // set face conn.
                pf->setPoint( 2, mesh.point( p2 ) ); // set face conn.
                pf->setPoint( 3, mesh.point( p3 ) ); // set face conn.
                pf->setPoint( 4, mesh.point( p4 ) ); // set face conn.
664
            }
665
666
667
668
669
670
            done++;
        }

        if ( line.find( "Edges" ) != std::string::npos )
        {
            nextIntINRIAMeshField( line.substr( line.find_last_of( "a" ) + 1 ), mystream );
671
            oStr << "Reading Bedges " << std::endl;
672
673
674
675
676
677
678
679
            for ( i = 0;i < nBEd;i++ )
            {
                mystream >> p1 >> p2 >> ibc;
                pe = &mesh.addEdge( true ); // Only boundary edges.
                pe->setMarker( EntityFlag( ibc ) );
                pe->setPoint( 1, mesh.point( p1 ) ); // set edge conn.
                pe->setPoint( 2, mesh.point( p2 ) ); // set edge conn.
            }
680
            oStr << "Boundary Edges Read " << std::endl;
681
682
683
684
685
686
            done++;
        }
        if ( line.find( "Tetrahedra" ) != std::string::npos )
        {
            count = 0;
            nextIntINRIAMeshField( line.substr( line.find_last_of( "a" ) + 1 ), mystream );
687
            oStr << "Reading Volumes " << std::endl;
688
689
690
691
            for ( i = 0; i < nVo; i++ )
            {
                mystream >> p1 >> p2 >> p3 >> p4 >> ibc;
                pv = &mesh.addVolume();
simone's avatar
simone committed
692
693
                pv->setId     ( i + 1 );
                pv->setLocalId( i + 1);
694
695
696
697
698
                pv->setPoint( 1, mesh.point( p1 ) );
                pv->setPoint( 2, mesh.point( p2 ) );
                pv->setPoint( 3, mesh.point( p3 ) );
                pv->setPoint( 4, mesh.point( p4 ) );
                pv->setMarker( EntityFlag( ibc ) );
699
700
//                mesh.localToGlobalElem().insert(std::make_pair(i+1, i+1));
//                mesh.globalToLocalElem().insert(std::make_pair(i+1, i+1));
701
702
                count++;
            }
703
704
            oStr << "size of the volume storage is "
                 << sizeof(VolumeType)*count/1024./1024. << " Mo." << std::endl;
705
            oStr << count << " Volume elements Read" << std::endl;
706

707
708
709
710
711
712
            done++;
        }
        if ( line.find( "Hexahedra" ) != std::string::npos )
        {
            count = 0;
            nextIntINRIAMeshField( line.substr( line.find_last_of( "a" ) + 1 ), mystream );
713
            oStr << "Reading Volumes " << std::endl;
714
715
716
717
            for ( i = 0; i < nVo; i++ )
            {
                mystream >> p1 >> p2 >> p3 >> p4 >> p5 >> p6 >> p7 >> p8 >> ibc;
                pv = &mesh.addVolume();
simone's avatar
simone committed
718
719
720
                pv->setId     ( i + 1 );
                pv->setLocalId( i + 1);
//                pv->id() = i + 1;
721
722
723
724
725
726
727
728
729
730
731
                pv->setPoint( 1, mesh.point( p1 ) );
                pv->setPoint( 2, mesh.point( p2 ) );
                pv->setPoint( 3, mesh.point( p3 ) );
                pv->setPoint( 4, mesh.point( p4 ) );
                pv->setPoint( 5, mesh.point( p5 ) );
                pv->setPoint( 6, mesh.point( p6 ) );
                pv->setPoint( 7, mesh.point( p7 ) );
                pv->setPoint( 8, mesh.point( p8 ) );
                pv->setMarker( EntityFlag( ibc ) );
                count++;
            }
732
            oStr << count << " Volume elements Read" << std::endl;
733
734
            done++;
        }
prudhomm's avatar
prudhomm committed
735

736
    }
737

738
739
    // Test mesh
    Switch sw;
prudhomm's avatar
prudhomm committed
740

741
    if ( !checkMesh3D( mesh, sw, true, verbose, oStr, std::cerr, oStr ) )
simone's avatar
simone committed
742
         abort();
743
    // if(!checkMesh3D(mesh, sw, true,true, oStr,oStr,oStr)) abort();//verbose version
prudhomm's avatar
prudhomm committed
744

745
    // This part is to build a P2 mesh from a P1 geometry
prudhomm's avatar
prudhomm committed
746

747
748
749
    if ( shape == TETRA && GeoShape::numPoints > 4 )
        p1top2( mesh );
    mystream.close();
750

751
    Real vols[ 3 ];
752
753
    getVolumeFromFaces( mesh, vols, oStr );
    oStr << "   VOLUME ENCLOSED BY THE MESH COMPUTED BY INTEGRATION ON" <<
754
        " BOUNDARY FACES" << std::endl;
755
    oStr << "INT(X)     INT(Y)      INT(Z) <- they should be equal and equal to" << std::endl <<
756
        "                                 the voulume enclosed by the mesh " << std::endl;
757
    oStr << vols[ 0 ] << " " << vols[ 1 ] << " " << vols[ 2 ] << std::endl;
prudhomm's avatar
prudhomm committed
758

759
    oStr << "   BOUNDARY FACES ARE DEFINING A CLOSED SURFACE IF " << testClosedDomain( mesh, oStr ) << std::endl <<
760
        " IS (ALMOST) ZERO" << std::endl;
prudhomm's avatar
prudhomm committed
761

762
    return done == 4 ;
prudhomm's avatar
prudhomm committed
763

764
}
simone's avatar
simone committed
765
766


prudhomm's avatar
prudhomm committed
767
//
768
769
// GMSH
//
770

prudhomm's avatar
prudhomm committed
771
/**
772
773
774
775
   read a gmsh mesh (3D) file and store it in a RegionMesh3D
   @param mesh mesh data structure to fill in
   @param filename name of the gmsh mesh file  to read
   @param regionFlag identifier for the region
776
   @param verbose whether the function shall be verbose
777
   @return true if everything went fine, false otherwise
prudhomm's avatar
prudhomm committed
778
*/
779
template <typename GeoShape, typename MC>
780
bool
781
readGmshFile( RegionMesh3D<GeoShape, MC> & mesh,
782
              const std::string & filename,
783
784
              EntityFlag regionFlag,
              bool verbose=false )
785
{
786
    std::ifstream __is ( filename.c_str() );
simone's avatar
simone committed
787
    Debug() << "gmsh reading: "<< filename << "\n";
788

789
790
791
    char __buf[256];
    __is >> __buf;
    Debug() << "buf: "<< __buf << "\n";
simone's avatar
simone committed
792
    UInt __n;
793
794
    __is >> __n;
    Debug() << "number of nodes: " << __n;
795

796
797
    // Add Marker to list of Markers
    mesh.setMarker( regionFlag );
798

799
800
801

    std::vector<double> __x(3*__n);
    std::vector<bool> __isonboundary(__n);
simone's avatar
simone committed
802
    std::vector<UInt> __whichboundary(__n);
803
804
    Debug() << "reading "<< __n << " nodes\n";
    std::map<int,int> itoii;
simone's avatar
simone committed
805
    for( UInt __i = 0; __i < __n;++__i )
806
    {
simone's avatar
simone committed
807
        UInt __ni;
prudhomm's avatar
prudhomm committed
808
809
810
811
        __is >> __ni
             >> __x[3*__i]
             >> __x[3*__i+1]
             >> __x[3*__i+2];
812
813
814
815
816
817
818

        itoii[__ni-1] = __i;
    }
    __is >> __buf;
    Debug() << "buf: "<< __buf << "\n";
    __is >> __buf;
    Debug() << "buf: "<< __buf << "\n";
simone's avatar
simone committed
819
    UInt __nele;
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
    __is >> __nele;

    typename RegionMesh3D<GeoShape, MC>::EdgeType * pe = 0;
    typename RegionMesh3D<GeoShape, MC>::FaceType * pf = 0;
    typename RegionMesh3D<GeoShape, MC>::VolumeType * pv = 0;




    Debug() << "number of elements: " << __nele << "\n";
    std::vector<std::vector<int> > __e(__nele);
    std::vector<int> __et(__nele);
    std::vector<int> __etype( __nele );
    std::vector<int> __gt(16);
    __gt.assign( 16, 0 );

simone's avatar
simone committed
836
    for( UInt __i = 0; __i < __nele;++__i )
837
    {
prudhomm's avatar
prudhomm committed
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
        int __ne, __t, __tag, __np, __dummy;
        __is >> __ne
             >> __t
             >> __tag
             >> __dummy
             >> __np;


        ++__gt[ __t ];
        __etype[__i] = __t;
        __et[__i] = __tag;
        __e[__i].resize( __np );
        int __p = 0;
        while ( __p != __np )
        {
            __is >> __e[__i][__p];
854
855
856
            __e[__i][__p] = itoii[ __e[__i][__p]-1];
            __e[__i][__p] += 1;

prudhomm's avatar
prudhomm committed
857
858
            ++__p;
        }
859
860
861
862
863
    }
    std::for_each( __gt.begin(), __gt.end(),  std::cout << boost::lambda::_1 << " " );
    std::cout << "\n";

    // Euler formulas
simone's avatar
simone committed
864
865
866
    UInt n_volumes = __gt[4];
    UInt n_faces_boundary = __gt[2];
    UInt n_faces_total = 2*n_volumes+(n_faces_boundary/2);
867

simone's avatar
simone committed
868
    mesh.setMaxNumGlobalPoints( __n );
869
870
    // Only Boundary Edges (in a next version I will allow for different choices)
    mesh.setMaxNumEdges( __gt[1] );
simone's avatar
simone committed
871
872
873
    mesh.setNumEdges   ( __gt[1] ); // Here the REAL number of edges (all of them)
    mesh.setNumBEdges  ( __gt[1] );
    mesh.setMaxNumGlobalEdges( __gt[1] );
874
875
876
    Debug() << "number of edges= " << __gt[1] << "\n";
    // Only Boundary Faces
    mesh.setMaxNumFaces( n_faces_total );
simone's avatar
simone committed
877
878
879
880
881
882
    mesh.setNumFaces   ( n_faces_total ); // Here the REAL number of edges (all of them)
    //mesh.setMaxNumFaces( n_faces_boundary );
    //mesh.setNumFaces   ( n_faces_boundary ); // Here the REAL number of edges (all of them)
    mesh.setNumBFaces  ( n_faces_boundary );
    mesh.setMaxNumGlobalFaces( n_faces_total );
    Debug() << "number of faces= " << n_faces_boundary << "\n";
883
    mesh.setMaxNumVolumes( n_volumes, true );
simone's avatar
simone committed
884
    mesh.setMaxNumGlobalVolumes( n_volumes );
885
886
887
888
    Debug() << "number of volumes= " << n_volumes << "\n";

    __isonboundary.assign( __n, false );
    __whichboundary.assign( __n, 0 );
simone's avatar
simone committed
889
    for( UInt __i = 0; __i < __nele;++__i )
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
    {
        switch( __etype[__i] )
        {
            // triangular faces (linear)
            case 2:
            {
                __isonboundary[ __e[__i][0]-1 ] = true;
                __isonboundary[ __e[__i][1]-1 ] = true;
                __isonboundary[ __e[__i][2]-1 ] = true;

                __whichboundary[__e[__i][0]-1 ] = __et[__i];
                __whichboundary[__e[__i][1]-1 ] = __et[__i];
                __whichboundary[__e[__i][2]-1 ] = __et[__i];
            }
        }
    }
    // add the point to the mesh
    typename RegionMesh3D<GeoShape, MC>::PointType * pp = 0;

    mesh.setMaxNumPoints( __n, true );
simone's avatar
simone committed
910
    mesh.setNumVertices ( __n );
simone's avatar
simone committed
911
    mesh.setNumBVertices( std::count( __isonboundary.begin(), __isonboundary.end(), true ) );
simone's avatar
simone committed
912
    mesh.setNumBPoints  ( mesh.numBVertices() );
913
914
915
916
917
918

    Debug() << "number of points : " << mesh.numPoints() << "\n";
    Debug() << "number of boundary points : " << mesh.numBPoints() << "\n";
    Debug() << "number of vertices : " << mesh.numVertices() << "\n";
    Debug() << "number of boundary vertices : " << mesh.numBVertices() << "\n";

simone's avatar
simone committed
919
    for( UInt __i = 0; __i < __n;++__i )
920
921
922
    {
        pp = &mesh.addPoint( __isonboundary[ __i ] );
        pp->setMarker( __whichboundary[__i] );
simone's avatar
simone committed
923
924
        pp->setId     ( __i + 1 );
        pp->setLocalId( __i + 1 );
925
926
927
        pp->x() = __x[3*__i];
        pp->y() = __x[3*__i+1];
        pp->z() = __x[3*__i+2];
928
929
        mesh.localToGlobalNode().insert(std::make_pair(__i+1, __i+1));
        mesh.globalToLocalNode().insert(std::make_pair(__i+1, __i+1));
930
931
    }

simone's avatar
simone committed
932
    int nVo = 1;
933
    // add the element to the mesh
simone's avatar
simone committed
934
    for( UInt __i = 0; __i < __nele;++__i )
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
    {
        switch( __etype[__i] )
        {
            // segment(linear)
            case 1:
            {
                pe = &( mesh.addEdge( true ) );
                pe->setMarker( EntityFlag( __et[__i] ) );
                pe->setPoint( 1, mesh.point( __e[__i][0] ) );
                pe->setPoint( 2, mesh.point( __e[__i][1] ) );



            }
            break;
            // triangular faces (linear)
            case 2:
            {
                pf = &( mesh.addFace( true ) );
                pf->setMarker( EntityFlag( __et[__i] ) );
                pf->setPoint( 1, mesh.point( __e[__i][0] ) );
                pf->setPoint( 2, mesh.point( __e[__i][1] ) );
                pf->setPoint( 3, mesh.point( __e[__i][2] ) );

            }
            break;
            // quadrangular faces(linear)
            case 3:
            {
                pf = &( mesh.addFace( true ) );
                pf->setMarker( EntityFlag( __et[__i] ) );
                pf->setPoint( 1, mesh.point( __e[__i][0] ) );
                pf->setPoint( 2, mesh.point( __e[__i][1] ) );
                pf->setPoint( 3, mesh.point( __e[__i][2] ) );
                pf->setPoint( 4, mesh.point( __e[__i][3] ) );
            }
            break;
            // tetrahedrons(linear)
            case 4:
            {
                pv = &( mesh.addVolume() );
simone's avatar
simone committed
976
977
978
                pv->setId     ( nVo );
                pv->setLocalId( nVo++);
//                pv->id() = __i + 1;
979
980
981
982
983
984
985
986
987
988
989
990
                pv->setMarker( EntityFlag( __et[__i] ) );
                pv->setPoint( 1, mesh.point( __e[__i][0] ) );
                pv->setPoint( 2, mesh.point( __e[__i][1] ) );
                pv->setPoint( 3, mesh.point( __e[__i][2] ) );
                pv->setPoint( 4, mesh.point( __e[__i][3] ) );
            }
            break;
            // hexahedrons(linear)
            case 5:
            {
                pv = &( mesh.addVolume() );

simone's avatar
simone committed
991
992
993
                pv->setId     ( __i + 1 );
                pv->setLocalId( __i + 1 );
//                pv->id() = __i + 1;
994
995
996
997
998
999
1000
                pv->setMarker( EntityFlag( __et[__i] ) );
                pv->setPoint( 1, mesh.point( __e[__i][0] ) );
                pv->setPoint( 2, mesh.point( __e[__i][1] ) );
                pv->setPoint( 3, mesh.point( __e[__i][2] ) );
                pv->setPoint( 4, mesh.point( __e[__i][3] ) );
                pv->setPoint( 5, mesh.point( __e[__i][4] ) );
                pv->setPoint( 6, mesh.point( __e[__i][5] ) );
For faster browsing, not all history is shown. View entire blame