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) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Reader module for Fieldview9 to read OpenFOAM mesh and data.
  26
  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 OpenFOAM 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 (isDir(rootAndCase/"system") && isDir(rootAndCase/"constant"))
 136     if (isDir(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.empty())
 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 an OpenFOAM 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.empty())
 540     {
 541         caseName = rootDir.name();
 542         rootDir  = rootDir.path();
 543     }
 544
 545     Info<< "rootDir  : " << rootDir << endl
 546         << "caseName : " << caseName << endl
 547         << "setName  : " << setName << endl;
 548
 549     //
 550     // Get/reuse database and mesh
 551     //
 552
 553     bool caseChanged = db_.setRunTime(rootDir, caseName, setName);
 554
 555
 556     //
 557     // Select time
 558     //
 559
 560     instantList Times = db_.runTime().times();
 561
 562     // If topo case set database time and update mesh.
 563     if (hasTopoChange(Times))
 564     {
 565         if (!selectTime(Times, iret))
 566         {
 567             return;
 568         }
 569     }
 570     else if (caseChanged)
 571     {
 572         // Load mesh (if case changed) to make sure we have nPoints etc.
 573         db_.loadMesh();
 574     }
 575
 576
 577     //
 578     // Set output variables
 579     //
 580
 581     *num_grids = 1;
 582
 583     const fvMesh& mesh = db_.mesh();
 584
 585     label nTotPoints = mesh.nPoints() + db_.polys().size();
 586
 587     num_nodes[0] = nTotPoints;
 588
 589     Info<< "setting num_nodes:" << num_nodes[0] << endl;
 590
 591     Info<< "setting num_face_types:" << mesh.boundary().size() << endl;
 592
 593     *num_face_types = mesh.boundary().size();
 594
 595     if (*num_face_types > *max_face_types)
 596     {
 597         errorMsg("Too many patches. FieldView limit:" + name(*max_face_types));
 598
 599         *iret = 1;
 600
 601         return;
 602     }
 603
 604
 605     forAll(mesh.boundaryMesh(), patchI)
 606     {
 607         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 608
 609         strcpy(face_type_names[patchI], patch.name().c_str());
 610
 611         if (isA<wallPolyPatch>(patch))
 612         {
 613             wall_flags[patchI] = 1;
 614         }
 615         else
 616         {
 617             wall_flags[patchI] = 0;
 618         }
 619         Info<< "Patch " << patch.name() << " is wall:"
 620             <<  wall_flags[patchI] << endl;
 621     }
 622
 623     //- Find all volFields and add them to database
 624     createFieldNames(db_.runTime(), Times, setName);
 625
 626     *num_vars = db_.volScalarNames().size() + 3*db_.volVectorNames().size();
 627
 628     if (*num_vars > *max_vars)
 629     {
 630         errorMsg("Too many variables. FieldView limit:" + name(*max_vars));
 631
 632         *iret = 1;
 633
 634         return;
 635     }
 636
 637
 638     int nameI = 0;
 639
 640     forAll(db_.volScalarNames(), i)
 641     {
 642         const word& fieldName = db_.volScalarNames()[i];
 643
 644         const word& fvName = db_.getFvName(fieldName);
 645
 646         strcpy(var_names[nameI++], fvName.c_str());
 647     }
 648
 649     forAll(db_.volVectorNames(), i)
 650     {
 651         const word& fieldName = db_.volVectorNames()[i];
 652
 653         const word& fvName = db_.getFvName(fieldName);
 654
 655         strcpy(var_names[nameI++], (fvName + "x;" + fvName).c_str());
 656         strcpy(var_names[nameI++], (fvName + "y").c_str());
 657         strcpy(var_names[nameI++], (fvName + "z").c_str());
 658     }
 659
 660     *iret = 0;
 661 }
 662
 663
 664 /*
 665  * Callback for reading timestep. Taken from user/user_unstruct_combined.f
 666  */
 667 void user_read_one_grid_function
 668 (
 669     int* iunit,     /* in: fortran unit number */
 670     int* igrid,     /* in: grid number to read */
 671     int* nodecnt,   /* in: number of nodes to read */
 672     float xyz[],    /* out: coordinates of nodes: x1..xN y1..yN z1..zN */
 673     int* num_vars,  /* in: number of results per node */
 674     float vars[],   /* out: values per node */
 675     int* iret       /* out: return value */
 676 )
 677 {
 678     fprintf(stderr, "\n** user_read_one_grid_function\n");
 679
 680     if (*igrid != 1)
 681     {
 682         errorMsg("Illegal grid number " + Foam::name(*igrid));
 683
 684         *iret = 1;
 685
 686         return;
 687     }
 688
 689     // Get current time
 690     instantList Times = db_.runTime().times();
 691
 692     // Set database time and update mesh.
 693     // Note: this should not be nessecary here. We already have the correct
 694     // time set and mesh loaded. This is only nessecary because Fieldview
 695     // otherwise thinks the case is non-transient.
 696     if (!selectTime(Times, iret))
 697     {
 698         return;
 699     }
 700
 701
 702     const fvMesh& mesh = db_.mesh();
 703
 704     // With mesh now loaded check for change in number of points.
 705     label nTotPoints = mesh.nPoints() + db_.polys().size();
 706
 707     if (*nodecnt != nTotPoints)
 708     {
 709         errorMsg
 710         (
 711             "nodecnt differs from number of points in mesh.\nnodecnt:"
 712           + Foam::name(*nodecnt)
 713           + "  mesh:"
 714           + Foam::name(nTotPoints)
 715         );
 716
 717         *iret = 1;
 718
 719         return;
 720     }
 721
 722
 723     if
 724     (
 725         *num_vars
 726      != (db_.volScalarNames().size() + 3*db_.volVectorNames().size())
 727     )
 728     {
 729         errorMsg("Illegal number of variables " + name(*num_vars));
 730
 731         *iret = 1;
 732
 733         return;
 734     }
 735
 736     //
 737     // Set coordinates
 738     //
 739
 740     const pointField& points = mesh.points();
 741
 742     int xIndex = 0;
 743     int yIndex = xIndex + nTotPoints;
 744     int zIndex = yIndex + nTotPoints;
 745
 746     // Add mesh points first.
 747     forAll(points, pointI)
 748     {
 749         xyz[xIndex++] = points[pointI].x();
 750         xyz[yIndex++] = points[pointI].y();
 751         xyz[zIndex++] = points[pointI].z();
 752     }
 753
 754     // Add cell centres of polys
 755     const pointField& ctrs = mesh.cellCentres();
 756
 757     const labelList& polys = db_.polys();
 758
 759     forAll(polys, i)
 760     {
 761         label cellI = polys[i];
 762
 763         xyz[xIndex++] = ctrs[cellI].x();
 764         xyz[yIndex++] = ctrs[cellI].y();
 765         xyz[zIndex++] = ctrs[cellI].z();
 766     }
 767
 768
 769     //
 770     // Define elements by calling fv routines
 771     //
 772
 773     static const cellModel& tet = *(cellModeller::lookup("tet"));
 774     static const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
 775     static const cellModel& pyr = *(cellModeller::lookup("pyr"));
 776     static const cellModel& prism = *(cellModeller::lookup("prism"));
 777     static const cellModel& wedge = *(cellModeller::lookup("wedge"));
 778     static const cellModel& hex = *(cellModeller::lookup("hex"));
 779     //static const cellModel& splitHex = *(cellModeller::lookup("splitHex"));
 780
 781     int tetVerts[4];
 782     int pyrVerts[5];
 783     int prismVerts[6];
 784     int hexVerts[8];
 785
 786     int tetFaces[4];
 787     int pyrFaces[5];
 788     int prismFaces[5];
 789     int hexFaces[6];
 790
 791     const cellShapeList& cellShapes = mesh.cellShapes();
 792     const faceList& faces = mesh.faces();
 793     const cellList& cells = mesh.cells();
 794     const labelList& owner = mesh.faceOwner();
 795     label nPoints = mesh.nPoints();
 796
 797     // Fieldview face_types array with all faces marked as internal.
 798     labelList internalFaces(6, 0);
 799
 800     // Mark all cells next to boundary so we don't have to calculate
 801     // wall_types for internal cells and can just pass internalFaces.
 802     boolList wallCell(mesh.nCells(), false);
 803
 804     label nWallCells = 0;
 805
 806     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 807     {
 808         label cellI = owner[faceI];
 809
 810         if (!wallCell[cellI])
 811         {
 812             wallCell[cellI] = true;
 813
 814             nWallCells++;
 815         }
 816     }
 817
 818     label nPolys = 0;
 819
 820     forAll(cellShapes, cellI)
 821     {
 822         const cellShape& cellShape = cellShapes[cellI];
 823         const cellModel& cellModel = cellShape.model();
 824         const cell& cellFaces = cells[cellI];
 825
 826         int istat = 0;
 827
 828         if (cellModel == tet)
 829         {
 830             tetVerts[0] = cellShape[3] + 1;
 831             tetVerts[1] = cellShape[0] + 1;
 832             tetVerts[2] = cellShape[1] + 1;
 833             tetVerts[3] = cellShape[2] + 1;
 834
 835             if (wallCell[cellI])
 836             {
 837                 labelList faceTypes =
 838                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 839
 840                 tetFaces[0] = faceTypes[2];
 841                 tetFaces[1] = faceTypes[3];
 842                 tetFaces[2] = faceTypes[0];
 843                 tetFaces[3] = faceTypes[1];
 844
 845                 istat = create_tet(ITYPE, tetVerts, tetFaces);
 846             }
 847             else
 848             {
 849                 // All faces internal so use precalculated zero.
 850                 istat = create_tet(ITYPE, tetVerts, internalFaces.begin());
 851             }
 852         }
 853         else if (cellModel == tetWedge)
 854         {
 855             prismVerts[0] = cellShape[0] + 1;
 856             prismVerts[1] = cellShape[3] + 1;
 857             prismVerts[2] = cellShape[4] + 1;
 858             prismVerts[3] = cellShape[1] + 1;
 859             prismVerts[4] = cellShape[4] + 1;
 860             prismVerts[5] = cellShape[2] + 1;
 861
 862             if (wallCell[cellI])
 863             {
 864                 labelList faceTypes =
 865                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 866
 867                 prismFaces[0] = faceTypes[1];
 868                 prismFaces[1] = faceTypes[2];
 869                 prismFaces[2] = faceTypes[3];
 870                 prismFaces[3] = faceTypes[0];
 871                 prismFaces[4] = faceTypes[3];
 872
 873                 istat = create_prism(ITYPE, prismVerts, prismFaces);
 874             }
 875             else
 876             {
 877                 istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
 878             }
 879         }
 880         else if (cellModel == pyr)
 881         {
 882             pyrVerts[0] = cellShape[0] + 1;
 883             pyrVerts[1] = cellShape[1] + 1;
 884             pyrVerts[2] = cellShape[2] + 1;
 885             pyrVerts[3] = cellShape[3] + 1;
 886             pyrVerts[4] = cellShape[4] + 1;
 887
 888             if (wallCell[cellI])
 889             {
 890                 labelList faceTypes =
 891                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 892
 893                 pyrFaces[0] = faceTypes[0];
 894                 pyrFaces[1] = faceTypes[3];
 895                 pyrFaces[2] = faceTypes[2];
 896                 pyrFaces[3] = faceTypes[1];
 897                 pyrFaces[4] = faceTypes[4];
 898
 899                 istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
 900             }
 901             else
 902             {
 903                 istat = create_pyramid(ITYPE, pyrVerts, internalFaces.begin());
 904             }
 905         }
 906         else if (cellModel == prism)
 907         {
 908             prismVerts[0] = cellShape[0] + 1;
 909             prismVerts[1] = cellShape[3] + 1;
 910             prismVerts[2] = cellShape[4] + 1;
 911             prismVerts[3] = cellShape[1] + 1;
 912             prismVerts[4] = cellShape[5] + 1;
 913             prismVerts[5] = cellShape[2] + 1;
 914
 915             if (wallCell[cellI])
 916             {
 917                 labelList faceTypes =
 918                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 919
 920                 prismFaces[0] = faceTypes[4];
 921                 prismFaces[1] = faceTypes[2];
 922                 prismFaces[2] = faceTypes[3];
 923                 prismFaces[3] = faceTypes[0];
 924                 prismFaces[4] = faceTypes[1];
 925
 926                 istat = create_prism(ITYPE, prismVerts, prismFaces);
 927             }
 928             else
 929             {
 930                 istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
 931             }
 932         }
 933         else if (cellModel == wedge)
 934         {
 935             hexVerts[0] = cellShape[0] + 1;
 936             hexVerts[1] = cellShape[1] + 1;
 937             hexVerts[2] = cellShape[0] + 1;
 938             hexVerts[3] = cellShape[2] + 1;
 939             hexVerts[4] = cellShape[3] + 1;
 940             hexVerts[5] = cellShape[4] + 1;
 941             hexVerts[6] = cellShape[6] + 1;
 942             hexVerts[7] = cellShape[5] + 1;
 943
 944             if (wallCell[cellI])
 945             {
 946                 labelList faceTypes =
 947                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 948
 949                 hexFaces[0] = faceTypes[2];
 950                 hexFaces[1] = faceTypes[3];
 951                 hexFaces[2] = faceTypes[0];
 952                 hexFaces[3] = faceTypes[1];
 953                 hexFaces[4] = faceTypes[4];
 954                 hexFaces[5] = faceTypes[5];
 955
 956                 istat = create_hex(ITYPE, hexVerts, hexFaces);
 957             }
 958             else
 959             {
 960                 istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
 961             }
 962         }
 963         else if (cellModel == hex)
 964         {
 965             hexVerts[0] = cellShape[0] + 1;
 966             hexVerts[1] = cellShape[1] + 1;
 967             hexVerts[2] = cellShape[3] + 1;
 968             hexVerts[3] = cellShape[2] + 1;
 969             hexVerts[4] = cellShape[4] + 1;
 970             hexVerts[5] = cellShape[5] + 1;
 971             hexVerts[6] = cellShape[7] + 1;
 972             hexVerts[7] = cellShape[6] + 1;
 973
 974             if (wallCell[cellI])
 975             {
 976                 labelList faceTypes =
 977                     getFaceTypes(mesh, cellFaces, cellShape.faces());
 978
 979                 hexFaces[0] = faceTypes[0];
 980                 hexFaces[1] = faceTypes[1];
 981                 hexFaces[2] = faceTypes[4];
 982                 hexFaces[3] = faceTypes[5];
 983                 hexFaces[4] = faceTypes[2];
 984                 hexFaces[5] = faceTypes[3];
 985
 986                 istat = create_hex(ITYPE, hexVerts, hexFaces);
 987             }
 988             else
 989             {
 990                 istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
 991             }
 992         }
 993         else
 994         {
 995             forAll(cellFaces, cFaceI)
 996             {
 997                 label faceI = cellFaces[cFaceI];
 998
 999                 // Get Fieldview facetype (internal/on patch)
1000                 label fvFaceType = getFvType(mesh, faceI);
1001
1002                 const face& f = faces[faceI];
1003
1004                 // Indices into storage
1005                 label nQuads = 0;
1006                 label nTris = 0;
1007
1008                 // Storage for triangles and quads created by face
1009                 // decomposition (sized for worst case)
1010                 faceList quadFaces((f.size() - 2)/2);
1011                 faceList triFaces(f.size() - 2);
1012
1013                 f.trianglesQuads
1014                 (
1015                     points,
1016                     nTris,
1017                     nQuads,
1018                     triFaces,
1019                     quadFaces
1020                 );
1021
1022                 // Label of cell centre in fv point list.
1023                 label polyCentrePoint = nPoints + nPolys;
1024
1025                 for (label i=0; i<nTris; i++)
1026                 {
1027                     if (cellI == owner[faceI])
1028                     {
1029                         tetVerts[0] = triFaces[i][0] + 1;
1030                         tetVerts[1] = triFaces[i][1] + 1;
1031                         tetVerts[2] = triFaces[i][2] + 1;
1032                         tetVerts[3] = polyCentrePoint + 1;
1033                     }
1034                     else
1035                     {
1036                         tetVerts[0] = triFaces[i][2] + 1;
1037                         tetVerts[1] = triFaces[i][1] + 1;
1038                         tetVerts[2] = triFaces[i][0] + 1;
1039                         tetVerts[3] = polyCentrePoint + 1;
1040                     }
1041
1042                     if (wallCell[cellI])
1043                     {
1044                         // Outside face is one without polyCentrePoint
1045                         tetFaces[0] = fvFaceType;
1046                         tetFaces[1] = 0;
1047                         tetFaces[2] = 0;
1048                         tetFaces[3] = 0;
1049
1050                         istat = create_tet(ITYPE, tetVerts, tetFaces);
1051                     }
1052                     else
1053                     {
1054                         istat =
1055                             create_tet
1056                             (
1057                                 ITYPE,
1058                                 tetVerts,
1059                                 internalFaces.begin()
1060                             );
1061                     }
1062                 }
1063
1064                 for (label i=0; i<nQuads; i++)
1065                 {
1066                     if (cellI == owner[faceI])
1067                     {
1068                         pyrVerts[0] = quadFaces[i][3] + 1;
1069                         pyrVerts[1] = quadFaces[i][2] + 1;
1070                         pyrVerts[2] = quadFaces[i][1] + 1;
1071                         pyrVerts[3] = quadFaces[i][0] + 1;
1072                         pyrVerts[4] = polyCentrePoint + 1;
1073
1074                     }
1075                     else
1076                     {
1077                         pyrVerts[0] = quadFaces[i][0] + 1;
1078                         pyrVerts[1] = quadFaces[i][1] + 1;
1079                         pyrVerts[2] = quadFaces[i][2] + 1;
1080                         pyrVerts[3] = quadFaces[i][3] + 1;
1081                         pyrVerts[4] = polyCentrePoint + 1;
1082                     }
1083
1084                     if (wallCell[cellI])
1085                     {
1086                         // Outside face is one without polyCentrePoint
1087                         pyrFaces[0] = fvFaceType;
1088                         pyrFaces[1] = 0;
1089                         pyrFaces[2] = 0;
1090                         pyrFaces[3] = 0;
1091                         pyrFaces[4] = 0;
1092
1093                         istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
1094                     }
1095                     else
1096                     {
1097                         istat =
1098                             create_pyramid
1099                             (
1100                                 ITYPE,
1101                                 pyrVerts,
1102                                 internalFaces.begin()
1103                             );
1104                     }
1105                 }
1106
1107                 if (istat != 0)
1108                 {
1109                     errorMsg("Error during adding cell " + name(cellI));
1110
1111                     *iret = 1;
1112
1113                     return;
1114                 }
1115             }
1116
1117             nPolys++;
1118         }
1119
1120         if (istat != 0)
1121         {
1122             errorMsg("Error during adding cell " + name(cellI));
1123
1124             *iret = 1;
1125
1126             return;
1127         }
1128     }
1129
1130
1131
1132     //
1133     // Set fieldvalues
1134     //
1135
1136     pointMesh pMesh(mesh);
1137
1138     volPointInterpolation pInterp(mesh, pMesh);
1139
1140     int pointI = 0;
1141
1142     forAll(db_.volScalarNames(), i)
1143     {
1144         const word& fieldName = db_.volScalarNames()[i];
1145
1146         storeScalarField(pInterp, fieldName, vars, pointI);
1147     }
1148
1149
1150     forAll(db_.volVectorNames(), i)
1151     {
1152         const word& fieldName = db_.volVectorNames()[i];
1153
1154         storeVectorField(pInterp, fieldName, vars, pointI);
1155     }
1156
1157     // Return without failure.
1158     *iret = 0;
1159 }
1160
1161
1162 void register_data_readers()
1163 {
1164     /*
1165     **
1166     ** You should edit this file to "register" a user-defined data
1167     ** reader with FIELDVIEW, if the user functions being registered
1168     ** here are written in C.
1169     ** You should edit "ftn_register_data_readers.f" if the user functions
1170     ** being registered are written in Fortran.
1171     ** In either case, the user functions being registered may call other
1172     ** functions written in either language (C or Fortran); only the
1173     ** language of the "query" and "read" functions referenced here matters
1174     ** to FIELDVIEW.
1175     **
1176     ** The following shows a sample user-defined data reader being
1177     ** "registered" with FIELDVIEW.
1178     **
1179     ** The "extern void" declarations should match the names of the
1180     ** query and grid-reading functions you are providing. It is
1181     ** strongly suggested that all such names begin with "user" so
1182     ** as not to conflict with global names in FIELDVIEW.
1183     **
1184     ** You may call any combination of the data reader registration
1185     ** functions shown below ("register_two_file_reader" and/or
1186     ** "register_single_file_reader" and/or "register_single_unstruct_reader"
1187     ** and/or "register_double_unstruct_reader") as many times as you like,
1188     ** in order to create several different data readers. Each data reader
1189     ** should of course have different "query" and "read" functions, all of
1190     ** which should also appear in "extern void" declarations.
1191     **
1192     */
1193
1194     /*
1195     ** like this for combined unstructured grids & results in a single file
1196     */
1197     reg_single_unstruct_reader
1198     (
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 // ************************************************************************* //