applications/utilities/postProcessing/graphics/fieldview9Reader/fieldview9Reader.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Reader module for Fieldview9 to read Foam mesh&data.
  27     Creates new 'fvbin' type executable which needs to be installed in place
  28     of bin/fvbin.
  29
  30     Implements a reader for combined mesh&results on an unstructured mesh.
  31
  32     See Fieldview Release 9 Reference Manual and coding in user/ directory
  33     of the Fieldview release.
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "fvCFD.H"
  38 #include "IOobjectList.H"
  39 #include "GeometricField.H"
  40 #include "pointFields.H"
  41 #include "volPointInterpolation.H"
  42 #include "readerDatabase.H"
  43 #include "wallPolyPatch.H"
  44 #include "ListOps.H"
  45 #include "cellModeller.H"
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 extern "C"
  50 {
  51 // Define various Fieldview constants and prototypes
  52
  53 #include "fv_reader_tags.h"
  54
  55 static const int FVHEX     = 2;
  56 static const int FVPRISM   = 3;
  57 static const int FVPYRAMID = 4;
  58 static const int FVTET     = 1;
  59 static const int ITYPE     = 1;
  60
  61 unsigned int fv_encode_elem_header(int elem_type, int wall_info[]);
  62 void time_step_get_value(float*, int, int*, float*, int*);
  63 void fv_error_msg(const char*, const char*);
  64
  65 void reg_single_unstruct_reader
  66 (
  67     char *,
  68     void
  69     (
  70         char*, int*, int*, int*, int*, int*, int*,
  71         int[], int*, char[][80], int[], int*, char[][80], int*
  72     ),
  73     void
  74     (
  75         int*, int*, int*, float[], int*, float[], int*
  76     )
  77 );
  78
  79 int create_tet(const int, const int[8], const int[]);
  80 int create_pyramid(const int, const int[8], const int[]);
  81 int create_prism(const int, const int[8], const int[]);
  82 int create_hex(const int, const int[8], const int[]);
  83
  84 typedef unsigned char uChar;
  85 extern uChar create_bndry_face_buffered
  86 (
  87     int bndry_type,
  88     int num_verts,
  89     int verts[],
  90     int *normals_flags,
  91     int num_grid_nodes
  92 );
  93
  94 /*
  95  * just define empty readers here for ease of linking.
  96  * Comment out if you have doubly defined linking error on this symbol
  97  */
  98 void ftn_register_data_readers()
  99 {}
 100
 101 /*
 102  * just define empty readers here for ease of linking.
 103  * Comment out if you have doubly defined linking error on this symbol
 104  */
 105 void ftn_register_functions()
 106 {}
 107
 108 /*
 109  * just define empty readers here for ease of linking.
 110  * Comment out if you have doubly defined linking error on this symbol
 111  */
 112 void register_functions()
 113 {}
 114
 115
 116 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 117
 118
 119 //
 120 // Storage for all Foam state (mainly database & mesh)
 121 //
 122 static readerDatabase db_;
 123
 124
 125 // Write fv error message.
 126 static void errorMsg(const string& msg)
 127 {
 128     fv_error_msg("Foam Reader", msg.c_str());
 129 }
 130
 131
 132 // Simple check if directory is valid case directory.
 133 static bool validCase(const fileName& rootAndCase)
 134 {
 135     //if (dir(rootAndCase/"system") && dir(rootAndCase/"constant"))
 136     if (dir(rootAndCase/"constant"))
 137     {
 138         return true;
 139     }
 140     else
 141     {
 142         return false;
 143     }
 144 }
 145
 146
 147 // Check whether case has topo changes by looking back from last time
 148 // to first directory with polyMesh/cells.
 149 static bool hasTopoChange(const instantList& times)
 150 {
 151     label lastIndex = times.size()-1;
 152
 153     const Time& runTime = db_.runTime();
 154
 155     // Only set time; do not update mesh.
 156     runTime.setTime(times[lastIndex], lastIndex);
 157
 158     fileName facesInst(runTime.findInstance(polyMesh::meshSubDir, "faces"));
 159
 160     // See if cellInst is constant directory. Note extra .name() is for
 161     // parallel cases when runTime.constant is '../constant'
 162     if (facesInst != runTime.constant().name())
 163     {
 164         Info<< "Found cells file in " << facesInst << " so case has "
 165             << "topological changes" << endl;
 166
 167         return true;
 168     }
 169     else
 170     {
 171         Info<< "Found cells file in " << facesInst << " so case has "
 172             << "no topological changes" << endl;
 173
 174         return false;
 175     }
 176 }
 177
 178
 179 static bool selectTime(const instantList& times, int* iret)
 180 {
 181     List<float> fvTimes(2*times.size());
 182
 183     forAll(times, timeI)
 184     {
 185         fvTimes[2*timeI]   = float(timeI);
 186         fvTimes[2*timeI+1] = float(times[timeI].value());
 187     }
 188
 189     int istep;
 190
 191     float time;
 192
 193     *iret=0;
 194
 195     time_step_get_value(fvTimes.begin(), times.size(), &istep, &time, iret);
 196
 197     if (*iret == -5)
 198     {
 199         errorMsg("Out of memory.");
 200
 201         return false;
 202     }
 203     if (*iret == -15)
 204     {
 205         // Cancel action.
 206         return false;
 207     }
 208     if (*iret != 0)
 209     {
 210         errorMsg("Unspecified error.");
 211
 212         return false;
 213     }
 214     Info<< "Selected timeStep:" << istep << "  time:" << scalar(time) << endl;
 215
 216     // Set time and load mesh.
 217     db_.setTime(times[istep], istep);
 218
 219     return true;
 220 }
 221
 222
 223 // Gets (names of) all fields in all timesteps.
 224 static void createFieldNames
 225 (
 226     const Time& runTime,
 227     const instantList& Times,
 228     const word& setName
 229 )
 230 {
 231     // From foamToFieldView9/getFieldNames.H:
 232
 233     HashSet<word> volScalarHash;
 234     HashSet<word> volVectorHash;
 235     HashSet<word> surfScalarHash;
 236     HashSet<word> surfVectorHash;
 237
 238     if (setName.size() == 0)
 239     {
 240         forAll(Times, timeI)
 241         {
 242             const word& timeName = Times[timeI].name();
 243
 244             // Add all fields to hashtable
 245             IOobjectList objects(runTime, timeName);
 246
 247             wordList vsNames(objects.names(volScalarField::typeName));
 248
 249             forAll(vsNames, fieldI)
 250             {
 251                 volScalarHash.insert(vsNames[fieldI]);
 252             }
 253
 254             wordList vvNames(objects.names(volVectorField::typeName));
 255
 256             forAll(vvNames, fieldI)
 257             {
 258                 volVectorHash.insert(vvNames[fieldI]);
 259             }
 260         }
 261     }
 262     db_.setFieldNames(volScalarHash.toc(), volVectorHash.toc());
 263 }
 264
 265
 266 // Appends interpolated values of fieldName to vars array.
 267 static void storeScalarField
 268 (
 269     const volPointInterpolation& pInterp,
 270     const word& fieldName,
 271     float vars[],
 272     label& pointI
 273 )
 274 {
 275     label nPoints = db_.mesh().nPoints();
 276     label nTotPoints = nPoints + db_.polys().size();
 277
 278     // Check if present
 279     IOobject ioHeader
 280     (
 281         fieldName,
 282         db_.runTime().timeName(),
 283         db_.runTime(),
 284         IOobject::MUST_READ,
 285         IOobject::NO_WRITE
 286     );
 287
 288     if (ioHeader.headerOk())
 289     {
 290         Info<< "Storing " << nTotPoints << " of interpolated " << fieldName
 291             << endl;
 292
 293         volScalarField field(ioHeader, db_.mesh());
 294
 295         pointScalarField psf(pInterp.interpolate(field));
 296
 297         forAll(psf, i)
 298         {
 299             vars[pointI++] = float(psf[i]);
 300         }
 301
 302         const labelList& polys = db_.polys();
 303
 304         forAll(polys, i)
 305         {
 306             label cellI = polys[i];
 307
 308             vars[pointI++] = float(field[cellI]);
 309         }
 310     }
 311     else
 312     {
 313         Info<< "Storing " << nTotPoints << " of dummy values of " << fieldName
 314             << endl;
 315
 316         for(label i = 0; i < nPoints; i++)
 317         {
 318             vars[pointI++] = 0.0;
 319         }
 320
 321         const labelList& polys = db_.polys();
 322
 323         forAll(polys, i)
 324         {
 325             vars[pointI++] = 0.0;
 326         }
 327     }
 328 }
 329
 330
 331 // Appends interpolated values of fieldName to vars array.
 332 static void storeVectorField
 333 (
 334     const volPointInterpolation& pInterp,
 335     const word& fieldName,
 336     float vars[],
 337     label& pointI
 338 )
 339 {
 340     label nPoints = db_.mesh().nPoints();
 341
 342     label nTotPoints = nPoints + db_.polys().size();
 343
 344     // Check if present
 345     IOobject ioHeader
 346     (
 347         fieldName,
 348         db_.runTime().timeName(),
 349         db_.runTime(),
 350         IOobject::MUST_READ,
 351         IOobject::NO_WRITE
 352     );
 353
 354     if (ioHeader.headerOk())
 355     {
 356         Info<< "Storing " << nTotPoints << " of interpolated " << fieldName
 357             << endl;
 358
 359         volVectorField field(ioHeader, db_.mesh());
 360
 361         for (direction d = 0; d < vector::nComponents; d++)
 362         {
 363             tmp<volScalarField> tcomp = field.component(d);
 364             const volScalarField& comp = tcomp();
 365
 366             pointScalarField psf(pInterp.interpolate(comp));
 367
 368             forAll(psf, i)
 369             {
 370                 vars[pointI++] = float(psf[i]);
 371             }
 372
 373             const labelList& polys = db_.polys();
 374
 375             forAll(polys, i)
 376             {
 377                 label cellI = polys[i];
 378
 379                 vars[pointI++] = float(comp[cellI]);
 380             }
 381         }
 382     }
 383     else
 384     {
 385         Info<< "Storing " << nTotPoints << " of dummy values of " << fieldName
 386             << endl;
 387
 388         for (direction d = 0; d < vector::nComponents; d++)
 389         {
 390             for(label i = 0; i < nPoints; i++)
 391             {
 392                 vars[pointI++] = 0.0;
 393             }
 394
 395             const labelList& polys = db_.polys();
 396
 397             forAll(polys, i)
 398             {
 399                 vars[pointI++] = 0.0;
 400             }
 401         }
 402     }
 403 }
 404
 405
 406 // Returns Fieldview face_type of mesh face faceI.
 407 static label getFvType(const polyMesh& mesh, const label faceI)
 408 {
 409     return mesh.boundaryMesh().whichPatch(faceI) + 1;
 410 }
 411
 412
 413 // Returns Fieldview face_type of face f.
 414 static label getFaceType
 415 (
 416     const polyMesh& mesh,
 417     const labelList& faceLabels,
 418     const face& f
 419 )
 420 {
 421     // Search in subset faceLabels of faces for index of face f.
 422     const faceList& faces = mesh.faces();
 423
 424     forAll(faceLabels, i)
 425     {
 426         label faceI = faceLabels[i];
 427
 428         if (f == faces[faceI])
 429         {
 430             // Convert patch to Fieldview face_type.
 431             return getFvType(mesh, faceI);
 432         }
 433     }
 434
 435     FatalErrorIn("getFaceType")
 436         << "Cannot find face " << f << " in mesh face subset " << faceLabels
 437         << abort(FatalError);
 438
 439     return -1;
 440 }
 441
 442
 443 // Returns Fieldview face_types for set of faces
 444 static labelList getFaceTypes
 445 (
 446     const polyMesh& mesh,
 447     const labelList& cellFaces,
 448     const faceList& cellShapeFaces
 449 )
 450 {
 451     labelList faceLabels(cellShapeFaces.size());
 452
 453     forAll(cellShapeFaces, i)
 454     {
 455         faceLabels[i] = getFaceType(mesh, cellFaces, cellShapeFaces[i]);
 456     }
 457     return faceLabels;
 458 }
 459
 460
 461 /*
 462  * Callback for querying file contents. Taken from user/user_unstruct_combined.f
 463  */
 464 void user_query_file_function
 465 (
 466     /* input */
 467     char* fname,        /* filename */
 468     int* lenf,          /* length of fName */
 469     int* iunit,         /* fortran unit to use */
 470     int* max_grids,     /* maximum number of grids allowed */
 471     int* max_face_types,/* maximum number of face types allowed */
 472     int* max_vars,      /* maximum number of result variables allowed per */
 473                         /* grid point*/
 474
 475     /* output */
 476     int* num_grids,     /* number of grids that will be read from data file */
 477     int  num_nodes[],   /* (array of node counts for all grids) */
 478     int* num_face_types,        /* number of special face types */
 479     char face_type_names[][80], /* array of face-type names */
 480     int  wall_flags[],          /* array of flags for the "wall" behavior */
 481     int* num_vars,              /* number of result variables per grid point */
 482     char var_names[][80],       /* array of variable names */
 483     int* iret                   /* return value */
 484 )
 485 {
 486     fprintf(stderr, "\n** user_query_file_function\n");
 487
 488     string rootAndCaseString(fname, *lenf);
 489
 490     fileName rootAndCase(rootAndCaseString);
 491
 492     word setName("");
 493
 494     if (!validCase(rootAndCase))
 495     {
 496         setName = rootAndCase.name();
 497
 498         rootAndCase = rootAndCase.path();
 499
 500         word setDir = rootAndCase.name();
 501
 502         rootAndCase = rootAndCase.path();
 503
 504         word meshDir = rootAndCase.name();
 505
 506         rootAndCase = rootAndCase.path();
 507         rootAndCase = rootAndCase.path();
 508
 509         if
 510         (
 511             setDir == "sets"
 512          && meshDir == polyMesh::typeName
 513          && validCase(rootAndCase)
 514         )
 515         {
 516             // Valid set (hopefully - cannot check contents of setName yet).
 517         }
 518         else
 519         {
 520             errorMsg
 521             (
 522                 "Could not find system/ and constant/ directory in\n"
 523               + rootAndCase
 524               + "\nPlease select a Foam case directory."
 525             );
 526
 527             *iret = 1;
 528
 529             return;
 530         }
 531
 532     }
 533
 534     fileName rootDir(rootAndCase.path());
 535
 536     fileName caseName(rootAndCase.name());
 537
 538     // handle trailing '/'
 539     if (caseName.size() == 0)
 540     {
 541         caseName = rootDir.name();
 542         rootDir = rootDir.path();
 543     }
 544
 545
 546     Info<< "rootDir  : " << rootDir << endl
 547         << "caseName : " << caseName << endl
 548         << "setName  : " << setName << endl;
 549
 550     //
 551     // Get/reuse database and mesh
 552     //
 553
 554     bool caseChanged = db_.setRunTime(rootDir, caseName, setName);
 555
 556
 557     //
 558     // Select time
 559     //
 560
 561     instantList Times = db_.runTime().times();
 562
 563     // If topo case set database time and update mesh.
 564     if (hasTopoChange(Times))
 565     {
 566         if (!selectTime(Times, iret))
 567         {
 568             return;
 569         }
 570     }
 571     else if (caseChanged)
 572     {
 573         // Load mesh (if case changed) to make sure we have nPoints etc.
 574         db_.loadMesh();
 575     }
 576
 577
 578     //
 579     // Set output variables
 580     //
 581
 582     *num_grids = 1;
 583
 584     const fvMesh& mesh = db_.mesh();
 585
 586     label nTotPoints = mesh.nPoints() + db_.polys().size();
 587
 588     num_nodes[0] = nTotPoints;
 589
 590     Info<< "setting num_nodes:" << num_nodes[0] << endl;
 591
 592     Info<< "setting num_face_types:" << mesh.boundary().size() << endl;
 593
 594     *num_face_types = mesh.boundary().size();
 595
 596     if (*num_face_types > *max_face_types)
 597     {
 598         errorMsg("Too many patches. FieldView limit:" + name(*max_face_types));
 599
 600         *iret = 1;
 601
 602         return;
 603     }
 604
 605
 606     forAll(mesh.boundaryMesh(), patchI)
 607     {
 608         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 609
 610         strcpy(face_type_names[patchI], patch.name().c_str());
 611
 612         if (isType<wallPolyPatch>(patch))
 613         {
 614             wall_flags[patchI] = 1;
 615         }
 616         else
 617         {
 618             wall_flags[patchI] = 0;
 619         }
 620         Info<< "Patch " << patch.name() << " is wall:"
 621             <<  wall_flags[patchI] << endl;
 622     }
 623
 624     //- Find all volFields and add them to database
 625     createFieldNames(db_.runTime(), Times, setName);
 626
 627     *num_vars = db_.volScalarNames().size() + 3*db_.volVectorNames().size();
 628
 629     if (*num_vars > *max_vars)
 630     {
 631         errorMsg("Too many variables. FieldView limit:" + name(*max_vars));
 632
 633         *iret = 1;
 634
 635         return;
 636     }
 637
 638
 639     int nameI = 0;
 640
 641     forAll(db_.volScalarNames(), i)
 642     {
 643         const word& fieldName = db_.volScalarNames()[i];
 644
 645         const word& fvName = db_.getFvName(fieldName);
 646
 647         strcpy(var_names[nameI++], fvName.c_str());
 648     }
 649
 650     forAll(db_.volVectorNames(), i)
 651     {
 652         const word& fieldName = db_.volVectorNames()[i];
 653
 654         const word& fvName = db_.getFvName(fieldName);
 655
 656         strcpy(var_names[nameI++], (fvName + "x;" + fvName).c_str());
 657         strcpy(var_names[nameI++], (fvName + "y").c_str());
 658         strcpy(var_names[nameI++], (fvName + "z").c_str());
 659     }
 660
 661     *iret = 0;
 662 }
 663
 664
 665 /*
 666  * Callback for reading timestep. Taken from user/user_unstruct_combined.f
 667  */
 668 void user_read_one_grid_function
 669 (
 670     int* iunit,     /* in: fortran unit number */
 671     int* igrid,     /* in: grid number to read */
 672     int* nodecnt,   /* in: number of nodes to read */
 673     float xyz[],    /* out: coordinates of nodes: x1..xN y1..yN z1..zN */
 674     int* num_vars,  /* in: number of results per node */
 675     float vars[],   /* out: values per node */
 676     int* iret       /* out: return value */
 677 )
 678 {
 679     fprintf(stderr, "\n** user_read_one_grid_function\n");
 680
 681     if (*igrid != 1)
 682     {
 683         errorMsg("Illegal grid number " + Foam::name(*igrid));
 684
 685         *iret = 1;
 686
 687         return;
 688     }
 689
 690     // Get current time
 691     instantList Times = db_.runTime().times();
 692
 693     // Set database time and update mesh.
 694     // Note: this should not be nessecary here. We already have the correct
 695     // time set and mesh loaded. This is only nessecary because Fieldview
 696     // otherwise thinks the case is non-transient.
 697     if (!selectTime(Times, iret))
 698     {
 699         return;
 700     }
 701
 702
 703     const fvMesh& mesh = db_.mesh();
 704
 705     // With mesh now loaded check for change in number of points.
 706     label nTotPoints = mesh.nPoints() + db_.polys().size();
 707
 708     if (*nodecnt != nTotPoints)
 709     {
 710         errorMsg
 711         (
 712             "nodecnt differs from number of points in mesh.\nnodecnt:"
 713           + Foam::name(*nodecnt)
 714           + "  mesh:"
 715           + Foam::name(nTotPoints)
 716         );
 717
 718         *iret = 1;
 719
 720         return;
 721     }
 722
 723
 724     if
 725     (
 726         *num_vars
 727      != (db_.volScalarNames().size() + 3*db_.volVectorNames().size())
 728     )
 729     {
 730         errorMsg("Illegal number of variables " + name(*num_vars));
 731
 732         *iret = 1;
 733
 734         return;
 735     }
 736
 737     //
 738     // Set coordinates
 739     //
 740
 741     const pointField& points = mesh.points();
 742
 743     int xIndex = 0;
 744     int yIndex = xIndex + nTotPoints;
 745     int zIndex = yIndex + nTotPoints;
 746
 747     // Add mesh points first.
 748     forAll(points, pointI)
 749     {
 750         xyz[xIndex++] = points[pointI].x();
 751         xyz[yIndex++] = points[pointI].y();
 752         xyz[zIndex++] = points[pointI].z();
 753     }
 754
 755     // Add cell centres of polys
 756     const pointField& ctrs = mesh.cellCentres();
 757
 758     const labelList& polys = db_.polys();
 759
 760     forAll(polys, i)
 761     {
 762         label cellI = polys[i];
 763
 764         xyz[xIndex++] = ctrs[cellI].x();
 765         xyz[yIndex++] = ctrs[cellI].y();
 766         xyz[zIndex++] = ctrs[cellI].z();
 767     }
 768
 769
 770     //
 771     // Define elements by calling fv routines
 772     //
 773
 774     static const cellModel& tet = *(cellModeller::lookup("tet"));
 775     static const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
 776     static const cellModel& pyr = *(cellModeller::lookup("pyr"));
 777     static const cellModel& prism = *(cellModeller::lookup("prism"));
 778     static const cellModel& wedge = *(cellModeller::lookup("wedge"));
 779     static const cellModel& hex = *(cellModeller::lookup("hex"));
 780     //static const cellModel& splitHex = *(cellModeller::lookup("splitHex"));
 781
 782     int tetVerts[4];
 783     int pyrVerts[5];
 784     int prismVerts[6];
 785     int hexVerts[8];
 786
 787     int tetFaces[4];
 788     int pyrFaces[5];
 789     int prismFaces[5];
 790     int hexFaces[6];
 791
 792     const cellShapeList& cellShapes = mesh.cellShapes();
 793     const faceList& faces = mesh.faces();
 794     const cellList& cells = mesh.cells();
 795     const labelList& owner = mesh.faceOwner();
 796     label nPoints = mesh.nPoints();
 797
 798     // Fieldview face_types array with all faces marked as internal.
 799     labelList internalFaces(6, 0);
 800
 801     // Mark all cells next to boundary so we don't have to calculate
 802     // wall_types for internal cells and can just pass internalFaces.
 803     boolList wallCell(mesh.nCells(), false);
 804
 805     label nWallCells = 0;
 806
 807     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 808     {
 809         label cellI = owner[faceI];
 810
 811         if (!wallCell[cellI])
 812         {
 813             wallCell[cellI] = true;
 814
 815             nWallCells++;
 816         }
 817     }
 818
 819     label nPolys = 0;
 820
 821     forAll(cellShapes, cellI)
 822     {
 823         const cellShape& cellShape = cellShapes[cellI];
 824         const cellModel& cellModel = cellShape.model();
 825         const cell& cellFaces = cells[cellI];
 826
 827         int istat = 0;
 828
 829         if (cellModel == tet)
 830         {
 831             tetVerts[0] = cellShape[3] + 1;
 832             tetVerts[1] = cellShape[0] + 1;
 833             tetVerts[2] = cellShape[1] + 1;
 834             tetVerts[3] = cellShape[2] + 1;
 835
 836             if (wallCell[cellI])
 837             {
 838                 labelList faceTypes =
 839                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 840
 841                 tetFaces[0] = faceTypes[2];
 842                 tetFaces[1] = faceTypes[3];
 843                 tetFaces[2] = faceTypes[0];
 844                 tetFaces[3] = faceTypes[1];
 845
 846                 istat = create_tet(ITYPE, tetVerts, tetFaces);
 847             }
 848             else
 849             {
 850                 // All faces internal so use precalculated zero.
 851                 istat = create_tet(ITYPE, tetVerts, internalFaces.begin());
 852             }
 853         }
 854         else if (cellModel == tetWedge)
 855         {
 856             prismVerts[0] = cellShape[0] + 1;
 857             prismVerts[1] = cellShape[3] + 1;
 858             prismVerts[2] = cellShape[4] + 1;
 859             prismVerts[3] = cellShape[1] + 1;
 860             prismVerts[4] = cellShape[4] + 1;
 861             prismVerts[5] = cellShape[2] + 1;
 862
 863             if (wallCell[cellI])
 864             {
 865                 labelList faceTypes =
 866                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 867
 868                 prismFaces[0] = faceTypes[1];
 869                 prismFaces[1] = faceTypes[2];
 870                 prismFaces[2] = faceTypes[3];
 871                 prismFaces[3] = faceTypes[0];
 872                 prismFaces[4] = faceTypes[3];
 873
 874                 istat = create_prism(ITYPE, prismVerts, prismFaces);
 875             }
 876             else
 877             {
 878                 istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
 879             }
 880         }
 881         else if (cellModel == pyr)
 882         {
 883             pyrVerts[0] = cellShape[0] + 1;
 884             pyrVerts[1] = cellShape[1] + 1;
 885             pyrVerts[2] = cellShape[2] + 1;
 886             pyrVerts[3] = cellShape[3] + 1;
 887             pyrVerts[4] = cellShape[4] + 1;
 888
 889             if (wallCell[cellI])
 890             {
 891                 labelList faceTypes =
 892                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 893
 894                 pyrFaces[0] = faceTypes[0];
 895                 pyrFaces[1] = faceTypes[3];
 896                 pyrFaces[2] = faceTypes[2];
 897                 pyrFaces[3] = faceTypes[1];
 898                 pyrFaces[4] = faceTypes[4];
 899
 900                 istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
 901             }
 902             else
 903             {
 904                 istat = create_pyramid(ITYPE, pyrVerts, internalFaces.begin());
 905             }
 906         }
 907         else if (cellModel == prism)
 908         {
 909             prismVerts[0] = cellShape[0] + 1;
 910             prismVerts[1] = cellShape[3] + 1;
 911             prismVerts[2] = cellShape[4] + 1;
 912             prismVerts[3] = cellShape[1] + 1;
 913             prismVerts[4] = cellShape[5] + 1;
 914             prismVerts[5] = cellShape[2] + 1;
 915
 916             if (wallCell[cellI])
 917             {
 918                 labelList faceTypes =
 919                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 920
 921                 prismFaces[0] = faceTypes[4];
 922                 prismFaces[1] = faceTypes[2];
 923                 prismFaces[2] = faceTypes[3];
 924                 prismFaces[3] = faceTypes[0];
 925                 prismFaces[4] = faceTypes[1];
 926
 927                 istat = create_prism(ITYPE, prismVerts, prismFaces);
 928             }
 929             else
 930             {
 931                 istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
 932             }
 933         }
 934         else if (cellModel == wedge)
 935         {
 936             hexVerts[0] = cellShape[0] + 1;
 937             hexVerts[1] = cellShape[1] + 1;
 938             hexVerts[2] = cellShape[0] + 1;
 939             hexVerts[3] = cellShape[2] + 1;
 940             hexVerts[4] = cellShape[3] + 1;
 941             hexVerts[5] = cellShape[4] + 1;
 942             hexVerts[6] = cellShape[6] + 1;
 943             hexVerts[7] = cellShape[5] + 1;
 944
 945             if (wallCell[cellI])
 946             {
 947                 labelList faceTypes =
 948                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 949
 950                 hexFaces[0] = faceTypes[2];
 951                 hexFaces[1] = faceTypes[3];
 952                 hexFaces[2] = faceTypes[0];
 953                 hexFaces[3] = faceTypes[1];
 954                 hexFaces[4] = faceTypes[4];
 955                 hexFaces[5] = faceTypes[5];
 956
 957                 istat = create_hex(ITYPE, hexVerts, hexFaces);
 958             }
 959             else
 960             {
 961                 istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
 962             }
 963         }
 964         else if (cellModel == hex)
 965         {
 966             hexVerts[0] = cellShape[0] + 1;
 967             hexVerts[1] = cellShape[1] + 1;
 968             hexVerts[2] = cellShape[3] + 1;
 969             hexVerts[3] = cellShape[2] + 1;
 970             hexVerts[4] = cellShape[4] + 1;
 971             hexVerts[5] = cellShape[5] + 1;
 972             hexVerts[6] = cellShape[7] + 1;
 973             hexVerts[7] = cellShape[6] + 1;
 974
 975             if (wallCell[cellI])
 976             {
 977                 labelList faceTypes =
 978                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 979
 980                 hexFaces[0] = faceTypes[0];
 981                 hexFaces[1] = faceTypes[1];
 982                 hexFaces[2] = faceTypes[4];
 983                 hexFaces[3] = faceTypes[5];
 984                 hexFaces[4] = faceTypes[2];
 985                 hexFaces[5] = faceTypes[3];
 986
 987                 istat = create_hex(ITYPE, hexVerts, hexFaces);
 988             }
 989             else
 990             {
 991                 istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
 992             }
 993         }
 994         else
 995         {
 996             forAll(cellFaces, cFaceI)
 997             {
 998                 label faceI = cellFaces[cFaceI];
 999
1000                 // Get Fieldview facetype (internal/on patch)
1001                 label fvFaceType = getFvType(mesh, faceI);
1002
1003                 const face& f = faces[faceI];
1004
1005                 // Indices into storage
1006                 label nQuads = 0;
1007                 label nTris = 0;
1008
1009                 // Storage for triangles and quads created by face
1010                 // decomposition (sized for worst case)
1011                 faceList quadFaces((f.size() - 2)/2);
1012                 faceList triFaces(f.size() - 2);
1013
1014                 f.trianglesQuads
1015                 (
1016                     points,
1017                     nTris,
1018                     nQuads,
1019                     triFaces,
1020                     quadFaces
1021                 );
1022
1023                 // Label of cell centre in fv point list.
1024                 label polyCentrePoint = nPoints + nPolys;
1025
1026                 for (label i=0; i<nTris; i++)
1027                 {
1028                     if (cellI == owner[faceI])
1029                     {
1030                         tetVerts[0] = triFaces[i][0] + 1;
1031                         tetVerts[1] = triFaces[i][1] + 1;
1032                         tetVerts[2] = triFaces[i][2] + 1;
1033                         tetVerts[3] = polyCentrePoint + 1;
1034                     }
1035                     else
1036                     {
1037                         tetVerts[0] = triFaces[i][2] + 1;
1038                         tetVerts[1] = triFaces[i][1] + 1;
1039                         tetVerts[2] = triFaces[i][0] + 1;
1040                         tetVerts[3] = polyCentrePoint + 1;
1041                     }
1042
1043                     if (wallCell[cellI])
1044                     {
1045                         // Outside face is one without polyCentrePoint
1046                         tetFaces[0] = fvFaceType;
1047                         tetFaces[1] = 0;
1048                         tetFaces[2] = 0;
1049                         tetFaces[3] = 0;
1050
1051                         istat = create_tet(ITYPE, tetVerts, tetFaces);
1052                     }
1053                     else
1054                     {
1055                         istat =
1056                             create_tet
1057                             (
1058                                 ITYPE,
1059                                 tetVerts,
1060                                 internalFaces.begin()
1061                             );
1062                     }
1063                 }
1064
1065                 for (label i=0; i<nQuads; i++)
1066                 {
1067                     if (cellI == owner[faceI])
1068                     {
1069                         pyrVerts[0] = quadFaces[i][3] + 1;
1070                         pyrVerts[1] = quadFaces[i][2] + 1;
1071                         pyrVerts[2] = quadFaces[i][1] + 1;
1072                         pyrVerts[3] = quadFaces[i][0] + 1;
1073                         pyrVerts[4] = polyCentrePoint + 1;
1074
1075                     }
1076                     else
1077                     {
1078                         pyrVerts[0] = quadFaces[i][0] + 1;
1079                         pyrVerts[1] = quadFaces[i][1] + 1;
1080                         pyrVerts[2] = quadFaces[i][2] + 1;
1081                         pyrVerts[3] = quadFaces[i][3] + 1;
1082                         pyrVerts[4] = polyCentrePoint + 1;
1083                     }
1084
1085                     if (wallCell[cellI])
1086                     {
1087                         // Outside face is one without polyCentrePoint
1088                         pyrFaces[0] = fvFaceType;
1089                         pyrFaces[1] = 0;
1090                         pyrFaces[2] = 0;
1091                         pyrFaces[3] = 0;
1092                         pyrFaces[4] = 0;
1093
1094                         istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
1095                     }
1096                     else
1097                     {
1098                         istat =
1099                             create_pyramid
1100                             (
1101                                 ITYPE,
1102                                 pyrVerts,
1103                                 internalFaces.begin()
1104                             );
1105                     }
1106                 }
1107
1108                 if (istat != 0)
1109                 {
1110                     errorMsg("Error during adding cell " + name(cellI));
1111
1112                     *iret = 1;
1113
1114                     return;
1115                 }
1116             }
1117
1118             nPolys++;
1119         }
1120
1121         if (istat != 0)
1122         {
1123             errorMsg("Error during adding cell " + name(cellI));
1124
1125             *iret = 1;
1126
1127             return;
1128         }
1129     }
1130
1131
1132
1133     //
1134     // Set fieldvalues
1135     //
1136
1137     pointMesh pMesh(mesh);
1138
1139     volPointInterpolation pInterp(mesh, pMesh);
1140
1141     int pointI = 0;
1142
1143     forAll(db_.volScalarNames(), i)
1144     {
1145         const word& fieldName = db_.volScalarNames()[i];
1146
1147         storeScalarField(pInterp, fieldName, vars, pointI);
1148     }
1149
1150
1151     forAll(db_.volVectorNames(), i)
1152     {
1153         const word& fieldName = db_.volVectorNames()[i];
1154
1155         storeVectorField(pInterp, fieldName, vars, pointI);
1156     }
1157
1158     // Return without failure.
1159     *iret = 0;
1160 }
1161
1162
1163 void register_data_readers()
1164 {
1165     /*
1166     **
1167     ** You should edit this file to "register" a user-defined data
1168     ** reader with FIELDVIEW, if the user functions being registered
1169     ** here are written in C.
1170     ** You should edit "ftn_register_data_readers.f" if the user functions
1171     ** being registered are written in Fortran.
1172     ** In either case, the user functions being registered may call other
1173     ** functions written in either language (C or Fortran); only the
1174     ** language of the "query" and "read" functions referenced here matters
1175     ** to FIELDVIEW.
1176     **
1177     ** The following shows a sample user-defined data reader being
1178     ** "registered" with FIELDVIEW.
1179     **
1180     ** The "extern void" declarations should match the names of the
1181     ** query and grid-reading functions you are providing. It is
1182     ** strongly suggested that all such names begin with "user" so
1183     ** as not to conflict with global names in FIELDVIEW.
1184     **
1185     ** You may call any combination of the data reader registration
1186     ** functions shown below ("register_two_file_reader" and/or
1187     ** "register_single_file_reader" and/or "register_single_unstruct_reader"
1188     ** and/or "register_double_unstruct_reader") as many times as you like,
1189     ** in order to create several different data readers. Each data reader
1190     ** should of course have different "query" and "read" functions, all of
1191     ** which should also appear in "extern void" declarations.
1192     **
1193     */
1194
1195     /*
1196     ** like this for combined unstructured grids & results in a single file
1197     */
1198     reg_single_unstruct_reader (
1199         "Foam Reader",                 /* title you want for data reader */
1200         user_query_file_function,      /* whatever you called this */
1201         user_read_one_grid_function    /* whatever you called this */
1202         );
1203 }
1204
1205
1206
1207
1208 }
1209
1210
1211 // ************************************************************************* //