applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Reads CCM files as written by Prostar/ccm using ccm 2.6 (not 2.4)
  26
  27     - does polyhedral mesh
  28     - does not handle 'interfaces' (couples)
  29     - does not handle cyclics. Use createPatch to recreate these
  30     - does not do data
  31     - does patch names only if they are in the problem description
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "ListOps.H"
  36 #include "argList.H"
  37 #include "foamTime.H"
  38 #include "fvMesh.H"
  39 #include "volFields.H"
  40 #include "emptyPolyPatch.H"
  41 #include "symmetryPolyPatch.H"
  42 #include "wallPolyPatch.H"
  43 #include "SortableList.H"
  44 #include "cellSet.H"
  45
  46 #include <ccmio.h>
  47 #include <vector>
  48
  49 using namespace Foam;
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53 static char const kDefaultState[] = "default";
  54 static int const kVertOffset = 2;
  55
  56
  57 // Determine upper-triangular order for internal faces:
  58 labelList getInternalFaceOrder
  59 (
  60     const cellList& cells,
  61     const labelList& owner,
  62     const labelList& neighbour
  63 )
  64 {
  65     labelList oldToNew(owner.size(), -1);
  66
  67     // First unassigned face
  68     label newFaceI = 0;
  69
  70     forAll(cells, cellI)
  71     {
  72         const labelList& cFaces = cells[cellI];
  73
  74         SortableList<label> nbr(cFaces.size());
  75
  76         forAll(cFaces, i)
  77         {
  78             label faceI = cFaces[i];
  79
  80             label nbrCellI = neighbour[faceI];
  81
  82             if (nbrCellI != -1)
  83             {
  84                 // Internal face. Get cell on other side.
  85                 if (nbrCellI == cellI)
  86                 {
  87                     nbrCellI = owner[faceI];
  88                 }
  89
  90                 if (cellI < nbrCellI)
  91                 {
  92                     // CellI is master
  93                     nbr[i] = nbrCellI;
  94                 }
  95                 else
  96                 {
  97                     // nbrCell is master. Let it handle this face.
  98                     nbr[i] = -1;
  99                 }
 100             }
 101             else
 102             {
 103                 // External face. Do later.
 104                 nbr[i] = -1;
 105             }
 106         }
 107
 108         nbr.sort();
 109
 110         forAll(nbr, i)
 111         {
 112             if (nbr[i] != -1)
 113             {
 114                 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
 115             }
 116         }
 117     }
 118
 119     // Keep boundary faces in same order.
 120     for (label faceI = newFaceI; faceI < owner.size(); faceI++)
 121     {
 122         oldToNew[faceI] = faceI;
 123     }
 124
 125     return oldToNew;
 126 }
 127
 128
 129 void storeCellInZone
 130 (
 131     const label cellI,
 132     const label cellType,
 133     Map<label>& typeToZone,
 134     List<DynamicList<label> >& zoneCells
 135 )
 136 {
 137     if (cellType >= 0)
 138     {
 139         Map<label>::iterator zoneFnd = typeToZone.find(cellType);
 140
 141         if (zoneFnd == typeToZone.end())
 142         {
 143             // New type. Allocate zone for it.
 144             zoneCells.setSize(zoneCells.size() + 1);
 145
 146             label zoneI = zoneCells.size()-1;
 147
 148             Info<< "Mapping type " << cellType << " to Foam cellZone "
 149                 << zoneI << endl;
 150             typeToZone.insert(cellType, zoneI);
 151
 152             zoneCells[zoneI].append(cellI);
 153         }
 154         else
 155         {
 156             // Existing zone for type
 157             zoneCells[zoneFnd()].append(cellI);
 158         }
 159     }
 160 }
 161
 162
 163 void CheckError(CCMIOError const &err, const Foam::string& str)
 164 {
 165     if (err != kCCMIONoErr)
 166     {
 167         FatalErrorIn("CheckError")
 168             << str << exit(FatalError);
 169     }
 170 }
 171
 172
 173 void ReadVertices
 174 (
 175     CCMIOError &err,
 176     CCMIOID &vertices,
 177     labelList& foamPointMap,
 178     pointField& foamPoints
 179 )
 180 {
 181     // Read the vertices.  This involves reading both the vertex data and
 182     // the map, which maps the index into the data array with the ID number.
 183     // As we process the vertices we need to be sure to scale them by the
 184     // appropriate scaling factor.  The offset is just to show you can read
 185     // any chunk.  Normally this would be in a for loop.
 186
 187     CCMIOSize nVertices;
 188     CCMIOEntitySize(&err, vertices, &nVertices, NULL);
 189
 190     List<int> mapData(nVertices, 0);
 191     List<float> verts(3*nVertices, 0);
 192
 193     int offset = 0;
 194     int offsetPlusSize = offset+nVertices;
 195
 196     int dims = 1;
 197     float scale;
 198     CCMIOID mapID;
 199     CCMIOReadVerticesf(&err, vertices, &dims, &scale, &mapID, verts.begin(),
 200         offset, offsetPlusSize);
 201     CCMIOReadMap(&err, mapID, mapData.begin(), offset, offsetPlusSize);
 202
 203     //CCMIOSize size;
 204     //CCMIOEntityDescription(&err, vertices, &size, NULL);
 205     //char *desc = new char[size + 1];
 206     //CCMIOEntityDescription(&err, vertices, NULL, desc);
 207     //Pout<< "label: '" << desc << "'" << endl;
 208     //delete [] desc;
 209
 210     // Convert to foamPoints
 211     foamPoints.setSize(nVertices);
 212     foamPoints = vector::zero;
 213     foamPointMap.setSize(nVertices);
 214
 215     forAll(foamPointMap, i)
 216     {
 217         foamPointMap[i] = mapData[i];
 218         for (direction cmpt = 0; cmpt < dims; cmpt++)
 219         {
 220             foamPoints[i][cmpt] = verts[dims*i + cmpt]*scale;
 221         }
 222     }
 223 }
 224
 225
 226 void ReadProblem
 227 (
 228     CCMIOError &err,
 229     CCMIOID& problem,
 230
 231     const Map<label>& prostarToFoamPatch,
 232     Map<word>& foamCellTypeNames,
 233     wordList& foamPatchTypes,
 234     wordList& foamPatchNames
 235 )
 236 {
 237     // ... walk through each cell type and print it...
 238     CCMIOID next;
 239     int i = 0;
 240     while
 241     (
 242         CCMIONextEntity(NULL, problem, kCCMIOCellType, &i, &next)
 243      == kCCMIONoErr
 244     )
 245     {
 246         char *name;
 247         int size;
 248         int cellType;
 249
 250         // ... if it has a material type.  (Note that we do not pass in
 251         // an array to get the name because we do not know how long the
 252         // string is yet.  Many parameters to CCMIO functions that
 253         // return
 254         // data can be NULL if that data is not needed.)
 255         if
 256         (
 257             CCMIOReadOptstr(NULL, next, "MaterialType", &size, NULL)
 258          == kCCMIONoErr
 259         )
 260         {
 261             name = new char[size + 1];
 262             CCMIOReadOptstr(&err, next, "MaterialType", &size, name);
 263             CCMIOGetEntityIndex(&err, next, &cellType);
 264
 265             foamCellTypeNames.insert(cellType, name);
 266             Pout<< "Celltype:" << cellType << " name:" << name << endl;
 267
 268             delete [] name;
 269         }
 270     }
 271
 272     // ... walk through each region description and print it...
 273
 274
 275     CCMIOID boundary;
 276     label regionI = 0;
 277     int k = 0;
 278     while
 279     (
 280         CCMIONextEntity(NULL, problem, kCCMIOBoundaryRegion, &k, &boundary)
 281      == kCCMIONoErr
 282     )
 283     {
 284         // Index of foam patch
 285         label foamPatchI = -1;
 286
 287         // Read prostar id
 288
 289         int prostarI = -1;
 290         if
 291         (
 292             CCMIOReadOpti(NULL, boundary, "ProstarRegionNumber", &prostarI)
 293          == kCCMIONoErr
 294         )
 295         {
 296             Pout<< "For region:" << regionI
 297                 << " found ProstarRegionNumber:" << prostarI << endl;
 298         }
 299         else
 300         {
 301             prostarI = regionI;
 302
 303             Pout<< "For region:" << regionI
 304                 << "did not find ProstarRegionNumber entry. Assuming "
 305                 << prostarI << endl;
 306         }
 307
 308
 309         if (prostarToFoamPatch.found(prostarI))
 310         {
 311             foamPatchI = prostarToFoamPatch[prostarI];
 312
 313             // Read boundary type
 314
 315             int size;
 316             if
 317             (
 318                 CCMIOReadOptstr(NULL, boundary, "BoundaryType", &size, NULL)
 319              == kCCMIONoErr
 320             )
 321             {
 322                 char* s = new char[size + 1];
 323                 CCMIOReadOptstr(NULL, boundary, "BoundaryType", &size, s);
 324                 s[size] = '\0';
 325                 foamPatchTypes[foamPatchI] = string::validate<word>(string(s));
 326                 delete [] s;
 327             }
 328
 329
 330             //foamPatchMap.append(prostarI);
 331
 332
 333             // Read boundary name:
 334             // - from BoundaryName field (Prostar)
 335             // - from 'Label' field (ccm+?)
 336             // - make up one from prostar id.
 337
 338             if
 339             (
 340                 CCMIOReadOptstr(NULL, boundary, "BoundaryName", &size, NULL)
 341              == kCCMIONoErr
 342             )
 343             {
 344                 char* name = new char[size + 1];
 345                 CCMIOReadOptstr(NULL, boundary, "BoundaryName", &size, name);
 346                 name[size] = '\0';
 347                 foamPatchNames[foamPatchI] =
 348                     string::validate<word>(string(name));
 349
 350                 delete [] name;
 351             }
 352             else if
 353             (
 354                 CCMIOReadOptstr(NULL, boundary, "Label", &size, NULL)
 355              == kCCMIONoErr
 356             )
 357             {
 358                 char* name = new char[size + 1];
 359                 CCMIOReadOptstr(NULL, boundary, "Label", &size, name);
 360                 name[size] = '\0';
 361                 foamPatchNames[foamPatchI] =
 362                     string::validate<word>(string(name));
 363                 delete [] name;
 364             }
 365             else
 366             {
 367                 foamPatchNames[foamPatchI] =
 368                     foamPatchTypes[foamPatchI]
 369                   + Foam::name(foamPatchI);
 370                 Pout<< "Made up name:" << foamPatchNames[foamPatchI]
 371                     << endl;
 372             }
 373
 374             Pout<< "Read patch:" << foamPatchI
 375                 << " name:" << foamPatchNames[foamPatchI]
 376                 << " foamPatchTypes:" << foamPatchTypes[foamPatchI]
 377                 << endl;
 378         }
 379
 380         regionI++;
 381     }
 382 }
 383
 384
 385 void ReadCells
 386 (
 387     CCMIOError &err,
 388     CCMIOID& topology,
 389     labelList& foamCellMap,
 390     labelList& foamCellType,
 391     Map<label>& prostarToFoamPatch,
 392     DynamicList<label>& foamPatchSizes,
 393     DynamicList<label>& foamPatchStarts,
 394     labelList& foamFaceMap,
 395     labelList& foamOwner,
 396     labelList& foamNeighbour,
 397     faceList& foamFaces
 398 )
 399 {
 400     // Read the cells.
 401     // ~~~~~~~~~~~~~~~
 402
 403     //  Store cell IDs so that we know what cells are in
 404     // this post data.
 405     CCMIOID id;
 406     CCMIOGetEntity(&err, topology, kCCMIOCells, 0, &id);
 407     CCMIOSize nCells;
 408     CCMIOEntitySize(&err, id, &nCells, NULL);
 409
 410     std::vector<int> mapData(nCells);
 411     std::vector<int> cellType(nCells);
 412
 413     CCMIOID mapID;
 414     CCMIOReadCells(&err, id, &mapID, &cellType[0], 0, nCells);
 415     CCMIOReadMap(&err, mapID, &mapData[0], 0, nCells);
 416     CheckError(err, "Error reading cells");
 417
 418     foamCellMap.setSize(nCells);
 419     foamCellType.setSize(nCells);
 420     forAll(foamCellMap, i)
 421     {
 422         foamCellMap[i] = mapData[i];
 423         foamCellType[i] = cellType[i];
 424     }
 425
 426
 427     // Read the internal faces.
 428     // ~~~~~~~~~~~~~~~~~~~~~~~~
 429
 430     CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
 431     CCMIOSize nInternalFaces;
 432     CCMIOEntitySize(&err, id, &nInternalFaces, NULL);
 433     Pout<< "nInternalFaces:" << nInternalFaces << endl;
 434
 435     // Determine patch sizes before reading internal faces
 436     label foamNFaces = nInternalFaces;
 437     int index = 0;
 438     while
 439     (
 440         CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index, &id)
 441      == kCCMIONoErr
 442     )
 443     {
 444         CCMIOSize size;
 445         CCMIOEntitySize(&err, id, &size, NULL);
 446
 447         Pout<< "Read kCCMIOBoundaryFaces entry with " << size
 448             << " faces." << endl;
 449
 450         foamPatchStarts.append(foamNFaces);
 451         foamPatchSizes.append(size);
 452         foamNFaces += size;
 453     }
 454     foamPatchStarts.shrink();
 455     foamPatchSizes.shrink();
 456
 457     Pout<< "patchSizes:" << foamPatchSizes << endl;
 458     Pout<< "patchStarts:" << foamPatchStarts << endl;
 459     Pout<< "nFaces:" << foamNFaces << endl;
 460
 461     mapData.resize(nInternalFaces);
 462     CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
 463     CCMIOSize size;
 464     CCMIOReadFaces(&err, id, kCCMIOInternalFaces, NULL, &size, NULL,
 465         kCCMIOStart, kCCMIOEnd);
 466     std::vector<int> faces(size);
 467     CCMIOReadFaces(&err, id, kCCMIOInternalFaces, &mapID, NULL, &faces[0],
 468         kCCMIOStart, kCCMIOEnd);
 469     std::vector<int> faceCells(2*nInternalFaces);
 470     CCMIOReadFaceCells(&err, id, kCCMIOInternalFaces, &faceCells[0],
 471         kCCMIOStart, kCCMIOEnd);
 472     CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
 473     CheckError(err, "Error reading internal faces");
 474
 475     // Copy into Foam lists
 476     foamFaceMap.setSize(foamNFaces);
 477     foamFaces.setSize(foamNFaces);
 478     foamOwner.setSize(foamNFaces);
 479     foamNeighbour.setSize(foamNFaces);
 480
 481     unsigned int pos = 0;
 482
 483     for (unsigned int faceI = 0; faceI < nInternalFaces; faceI++)
 484     {
 485         foamFaceMap[faceI] = mapData[faceI];
 486         foamOwner[faceI] = faceCells[2*faceI];
 487         foamNeighbour[faceI] = faceCells[2*faceI+1];
 488         face& f = foamFaces[faceI];
 489
 490         f.setSize(faces[pos++]);
 491         forAll(f, fp)
 492         {
 493             f[fp] = faces[pos++];
 494         }
 495     }
 496
 497     // Read the boundary faces.
 498     // ~~~~~~~~~~~~~~~~~~~~~~~~
 499
 500     index = 0;
 501     label regionI = 0;
 502     while
 503     (
 504         CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index, &id)
 505      == kCCMIONoErr
 506     )
 507     {
 508         CCMIOSize nFaces;
 509         CCMIOEntitySize(&err, id, &nFaces, NULL);
 510
 511         mapData.resize(nFaces);
 512         faceCells.resize(nFaces);
 513         CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, NULL, &size, NULL,
 514                    kCCMIOStart, kCCMIOEnd);
 515         faces.resize(size);
 516         CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, &mapID, NULL, &faces[0],
 517                    kCCMIOStart, kCCMIOEnd);
 518         CCMIOReadFaceCells(&err, id, kCCMIOBoundaryFaces, &faceCells[0],
 519                    kCCMIOStart, kCCMIOEnd);
 520         CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
 521         CheckError(err, "Error reading boundary faces");
 522
 523         // Read prostar id
 524         int prostarI;
 525         if
 526         (
 527             CCMIOReadOpti(NULL, id, "ProstarRegionNumber", &prostarI)
 528          == kCCMIONoErr
 529         )
 530         {
 531             Pout<< "For region:" << regionI
 532                 << " found ProstarRegionNumber:" << prostarI << endl;
 533         }
 534         else
 535         {
 536             prostarI = regionI;
 537
 538             Pout<< "For region:" << regionI
 539                 << " did not find ProstarRegionNumber entry. Assuming "
 540                 << prostarI << endl;
 541         }
 542         prostarToFoamPatch.insert(prostarI, regionI);
 543
 544
 545         Pout<< "region:" << regionI
 546             << " ProstarRegionNumber:" << prostarI
 547             << " foamPatchStart:"
 548             << foamPatchStarts[regionI]
 549             << " size:"
 550             << foamPatchSizes[regionI]
 551             << endl;
 552
 553         // Copy into Foam list.
 554         unsigned int pos = 0;
 555
 556         for (unsigned int i = 0; i < nFaces; i++)
 557         {
 558             label foamFaceI = foamPatchStarts[regionI] + i;
 559
 560             foamFaceMap[foamFaceI] = mapData[i];
 561             foamOwner[foamFaceI] = faceCells[i];
 562             foamNeighbour[foamFaceI] = -1;
 563             face& f = foamFaces[foamFaceI];
 564
 565             f.setSize(faces[pos++]);
 566             forAll(f, fp)
 567             {
 568                 f[fp] = faces[pos++];
 569             }
 570         }
 571         Pout<< endl;
 572
 573         regionI++;
 574     }
 575 }
 576
 577
 578 // Main program:
 579
 580 int main(int argc, char *argv[])
 581 {
 582     argList::noParallel();
 583     argList::validArgs.append("ccm26 file");
 584
 585 #   include "setRootCase.H"
 586 #   include "createTime.H"
 587
 588
 589     // Foam mesh data
 590     // ~~~~~~~~~~~~~~
 591
 592     // Coordinates
 593     pointField foamPoints;
 594     labelList foamPointMap;
 595
 596     // Cell type
 597     labelList foamCellType;
 598     labelList foamCellMap;
 599
 600     // Patching info
 601     Map<label> prostarToFoamPatch;
 602     DynamicList<label> foamPatchSizes;
 603     DynamicList<label> foamPatchStarts;
 604     // Face connectivity
 605     labelList foamFaceMap;
 606     labelList foamOwner;
 607     labelList foamNeighbour;
 608     faceList foamFaces;
 609
 610     // Celltypes (not the names of the zones; just the material type)
 611     // and patch type names
 612     Map<word> foamCellTypeNames;
 613     wordList foamPatchTypes;
 614     wordList foamPatchNames;
 615
 616     {
 617         fileName ccmFile(args.additionalArgs()[0]);
 618
 619         if (!isFile(ccmFile))
 620         {
 621             FatalErrorIn(args.executable())
 622                 << "Cannot read file " << ccmFile
 623                 << exit(FatalError);
 624         }
 625
 626         word ccmExt = ccmFile.ext();
 627
 628         if (ccmExt != "ccm" && ccmExt != "ccmg")
 629         {
 630             FatalErrorIn(args.executable())
 631                 << "Illegal extension " << ccmExt << " for file " << ccmFile
 632                 << nl << "Allowed extensions are '.ccm', '.ccmg'"
 633                 << exit(FatalError);
 634         }
 635
 636         // Open the file.  Because we did not initialize 'err' we need to pass
 637         // in NULL (which always means kCCMIONoErr) and then assign the return
 638         // value to 'err'.).
 639         CCMIOID root;
 640         CCMIOError err =
 641             CCMIOOpenFile(NULL, ccmFile.c_str(), kCCMIORead, &root);
 642
 643         // We are going to assume that we have a state with a known name.
 644         // We could instead use CCMIONextEntity() to walk through all the
 645         // states in the file and present the list to the user for selection.
 646         CCMIOID state;
 647         int stateI = 0;
 648         CCMIONextEntity(&err, root, kCCMIOState, &stateI, &state);
 649         CheckError(err, "Error opening state");
 650
 651         unsigned int size;
 652         CCMIOEntityDescription(&err, state, &size, NULL);
 653         char *desc = new char[size + 1];
 654         CCMIOEntityDescription(&err, state, NULL, desc);
 655         Pout<< "Reading state '" << kDefaultState << "' (" << desc << ")"
 656             << endl;
 657         delete [] desc;
 658
 659         // Find the first processor (i has previously been initialized to 0)
 660         //  and read the mesh and solution information.
 661         int i = 0;
 662         CCMIOID processor;
 663         CCMIONextEntity(&err, state, kCCMIOProcessor, &i, &processor);
 664         CCMIOID solution, vertices, topology;
 665         CCMIOReadProcessor
 666         (
 667             &err,
 668             processor,
 669             &vertices,
 670             &topology,
 671             NULL,
 672             &solution
 673         );
 674
 675         if (err != kCCMIONoErr)
 676         {
 677             // Maybe no solution;  try again
 678             err = kCCMIONoErr;
 679             CCMIOReadProcessor
 680             (
 681                 &err,
 682                 processor,
 683                 &vertices,
 684                 &topology,
 685                 NULL,
 686                 NULL
 687             );
 688             if (err != kCCMIONoErr)
 689             {
 690                 FatalErrorIn(args.executable())
 691                     << "Could not read the file."
 692                     << exit(FatalError);
 693             }
 694         }
 695
 696         ReadVertices(err, vertices, foamPointMap, foamPoints);
 697
 698         Pout<< "nPoints:" << foamPoints.size() << endl
 699             << "bounding box:" << boundBox(foamPoints) << endl
 700             << endl;
 701
 702         ReadCells
 703         (
 704             err,
 705             topology,
 706             foamCellMap,
 707             foamCellType,
 708             prostarToFoamPatch,
 709             foamPatchSizes,
 710             foamPatchStarts,
 711             foamFaceMap,
 712             foamOwner,
 713             foamNeighbour,
 714             foamFaces
 715         );
 716
 717         Pout<< "nCells:" << foamCellMap.size() << endl
 718             << "nFaces:" << foamOwner.size() << endl
 719             << "nPatches:" << foamPatchStarts.size() << endl
 720             << "nInternalFaces:" << foamPatchStarts[0] << endl
 721             << endl;
 722
 723         // Create some default patch names/types. These will be overwritten
 724         // by any problem desciption (if it is there)
 725         foamPatchTypes.setSize(foamPatchStarts.size());
 726         foamPatchNames.setSize(foamPatchStarts.size());
 727
 728         forAll(foamPatchNames, i)
 729         {
 730             foamPatchNames[i] = word("patch") + Foam::name(i);
 731             foamPatchTypes[i] = "patch";
 732         }
 733
 734         // Get problem description
 735
 736         CCMIOID problem;
 737         int problemI = 0;
 738         CCMIONextEntity
 739         (
 740             &err,
 741             root,
 742             kCCMIOProblemDescription,
 743             &problemI,
 744             &problem
 745         );
 746         CheckError(err, "Error stepping to first problem description");
 747
 748         if (CCMIOIsValidEntity(problem))   // if we have a problem description
 749         {
 750             ReadProblem
 751             (
 752                 err,
 753                 problem,
 754                 prostarToFoamPatch,
 755
 756                 foamCellTypeNames,
 757                 foamPatchTypes,
 758                 foamPatchNames
 759             );
 760         }
 761
 762
 763         CCMIOCloseFile(&err, vertices);
 764         CCMIOCloseFile(&err, topology);
 765         CCMIOCloseFile(&err, solution);
 766         CCMIOCloseFile(&err, root);
 767     }
 768
 769
 770     Pout<< "foamPatchNames:" << foamPatchNames << endl;
 771
 772
 773     Pout<< "foamOwner : min:" << min(foamOwner)
 774         << " max:" << max(foamOwner)
 775         << nl
 776         << "foamNeighbour : min:" << min(foamNeighbour)
 777         << " max:" << max(foamNeighbour)
 778         << nl
 779         << "foamCellType : min:" << min(foamCellType)
 780         << " max:" << max(foamCellType)
 781         << nl << endl;
 782
 783
 784
 785     // We now have extracted all info from CCMIO:
 786     // - coordinates (points)
 787     // - face to point addressing (faces)
 788     // - face to cell addressing (owner, neighbour)
 789     // - cell based data (cellId)
 790
 791
 792     // Renumber vertex labels to Foam point labels
 793     {
 794         label maxCCMPointI = max(foamPointMap);
 795         labelList toFoamPoints(invert(maxCCMPointI+1, foamPointMap));
 796
 797         forAll(foamFaces, faceI)
 798         {
 799             inplaceRenumber(toFoamPoints, foamFaces[faceI]);
 800         }
 801     }
 802
 803     // Renumber cell labels
 804     {
 805         label maxCCMCellI = max(foamCellMap);
 806         labelList toFoamCells(invert(maxCCMCellI+1, foamCellMap));
 807
 808         inplaceRenumber(toFoamCells, foamOwner);
 809         inplaceRenumber(toFoamCells, foamNeighbour);
 810     }
 811
 812
 813     //
 814     // Basic mesh info complete. Now convert to Foam convention:
 815     // - owner < neighbour
 816     // - face vertices such that normal points away from owner
 817     // - order faces: upper-triangular for internal faces; boundary faces after
 818     //   internal faces
 819     //
 820
 821     // Set owner/neighbour so owner < neighbour
 822     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 823
 824     forAll(foamNeighbour, faceI)
 825     {
 826         label nbr = foamNeighbour[faceI];
 827         label own = foamOwner[faceI];
 828
 829         if (nbr >= foamCellType.size() || own >= foamCellType.size())
 830         {
 831             FatalErrorIn(args.executable())
 832                 << "face:" << faceI
 833                 << " nbr:" << nbr
 834                 << " own:" << own
 835                 << " nCells:" << foamCellType.size()
 836                 << exit(FatalError);
 837         }
 838
 839         if (nbr >= 0)
 840         {
 841             if (nbr < own)
 842             {
 843                 foamOwner[faceI] = foamNeighbour[faceI];
 844                 foamNeighbour[faceI] = own;
 845                 foamFaces[faceI] = foamFaces[faceI].reverseFace();
 846             }
 847         }
 848
 849
 850         // And check the face
 851         const face& f = foamFaces[faceI];
 852
 853         forAll(f, fp)
 854         {
 855             if (f[fp] < 0 || f[fp] >= foamPoints.size())
 856             {
 857                 FatalErrorIn(args.executable()) << "face:" << faceI
 858                     << " f:" << f << abort(FatalError);
 859             }
 860         }
 861     }
 862
 863
 864     // Do upper-triangular ordering
 865     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 866
 867     {
 868         // Create cells (inverse of face-to-cell addressing)
 869         cellList foamCells;
 870         primitiveMesh::calcCells
 871         (
 872             foamCells,
 873             foamOwner,
 874             foamNeighbour,
 875             foamCellType.size()
 876         );
 877
 878         // Determine face order for upper-triangular ordering
 879         labelList oldToNew
 880         (
 881             getInternalFaceOrder
 882             (
 883                 foamCells,
 884                 foamOwner,
 885                 foamNeighbour
 886             )
 887         );
 888
 889         // Reorder faces accordingly
 890         forAll(foamCells, cellI)
 891         {
 892             inplaceRenumber(oldToNew, foamCells[cellI]);
 893         }
 894
 895         // Reorder faces.
 896         inplaceReorder(oldToNew, foamFaces);
 897         inplaceReorder(oldToNew, foamOwner);
 898         inplaceReorder(oldToNew, foamNeighbour);
 899     }
 900
 901
 902     // Construct fvMesh
 903     // ~~~~~~~~~~~~~~~~
 904
 905     fvMesh mesh
 906     (
 907         IOobject
 908         (
 909             fvMesh::defaultRegion,
 910             runTime.constant(),
 911             runTime
 912         ),
 913         xferMove<pointField>(foamPoints),
 914         xferMove<faceList>(foamFaces),
 915         xferCopy<labelList>(foamOwner),
 916         xferMove<labelList>(foamNeighbour)
 917     );
 918
 919     // Create patches. Use patch types to determine what Foam types to generate.
 920     List<polyPatch*> newPatches(foamPatchNames.size());
 921
 922     label meshFaceI = foamPatchStarts[0];
 923
 924     forAll(newPatches, patchI)
 925     {
 926         const word& patchName = foamPatchNames[patchI];
 927         const word& patchType = foamPatchTypes[patchI];
 928
 929         Pout<< "Patch:" << patchName << " start at:" << meshFaceI
 930             << " size:" << foamPatchSizes[patchI]
 931             << " end at:" << meshFaceI+foamPatchSizes[patchI]
 932             << endl;
 933
 934         if (patchType == "wall")
 935         {
 936             newPatches[patchI] =
 937                 new wallPolyPatch
 938                 (
 939                     patchName,
 940                     foamPatchSizes[patchI],
 941                     meshFaceI,
 942                     patchI,
 943                     mesh.boundaryMesh()
 944                 );
 945         }
 946         else if (patchType == "symmetryplane")
 947         {
 948             newPatches[patchI] =
 949                 new symmetryPolyPatch
 950                 (
 951                     patchName,
 952                     foamPatchSizes[patchI],
 953                     meshFaceI,
 954                     patchI,
 955                     mesh.boundaryMesh()
 956                 );
 957         }
 958         else if (patchType == "empty")
 959         {
 960             // Note: not ccm name, introduced by us above.
 961             newPatches[patchI] =
 962                 new emptyPolyPatch
 963                 (
 964                     patchName,
 965                     foamPatchSizes[patchI],
 966                     meshFaceI,
 967                     patchI,
 968                     mesh.boundaryMesh()
 969                 );
 970         }
 971         else
 972         {
 973             // All other ccm types become straight polyPatch:
 974             // 'inlet', 'outlet', 'pressured'.
 975             newPatches[patchI] =
 976                 new polyPatch
 977                 (
 978                     patchName,
 979                     foamPatchSizes[patchI],
 980                     meshFaceI,
 981                     patchI,
 982                     mesh.boundaryMesh()
 983                 );
 984         }
 985
 986         meshFaceI += foamPatchSizes[patchI];
 987     }
 988
 989     if (meshFaceI != foamOwner.size())
 990     {
 991         FatalErrorIn(args.executable())
 992             << "meshFaceI:" << meshFaceI
 993             << " nFaces:" << foamOwner.size()
 994             << abort(FatalError);
 995     }
 996     mesh.addFvPatches(newPatches);
 997
 998
 999
1000     // Construct sets and zones from types
1001     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1002
1003     label maxType = max(foamCellType);
1004     label minType = min(foamCellType);
1005
1006     if (maxType > minType)
1007     {
1008         // From foamCellType physical region to Foam cellZone
1009         Map<label> typeToZone;
1010         // Storage for cell zones.
1011         List<DynamicList<label> > zoneCells(0);
1012
1013         forAll(foamCellType, cellI)
1014         {
1015             storeCellInZone
1016             (
1017                 cellI,
1018                 foamCellType[cellI],
1019                 typeToZone,
1020                 zoneCells
1021             );
1022         }
1023
1024         mesh.cellZones().clear();
1025         mesh.cellZones().setSize(typeToZone.size());
1026
1027         label nValidCellZones = 0;
1028
1029         forAllConstIter(Map<label>, typeToZone, iter)
1030         {
1031             label type = iter.key();
1032             label zoneI = iter();
1033
1034             zoneCells[zoneI].shrink();
1035
1036             word zoneName = "cellZone_" + name(type);
1037
1038             Info<< "Writing zone " << type
1039                 << " size " << zoneCells[zoneI].size()
1040                 << " to cellZone "
1041                 << zoneName << " and cellSet " << zoneName
1042                 << endl;
1043
1044             cellSet cset(mesh, zoneName, labelHashSet(zoneCells[zoneI]));
1045             cset.write();
1046
1047             mesh.cellZones().set
1048             (
1049                 nValidCellZones,
1050                 new cellZone
1051                 (
1052                     zoneName,
1053                     zoneCells[zoneI],
1054                     nValidCellZones,
1055                     mesh.cellZones()
1056                 )
1057             );
1058             nValidCellZones++;
1059         }
1060         mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
1061     }
1062
1063
1064     Info<< "Writing mesh to " << mesh.objectRegistry::objectPath()
1065         << "..." << nl << endl;
1066
1067
1068     // Construct field with calculated bc to hold Star cell Id.
1069     volScalarField cellIdField
1070     (
1071         IOobject
1072         (
1073             "cellId",
1074             runTime.timeName(),
1075             mesh,
1076             IOobject::NO_READ,
1077             IOobject::AUTO_WRITE
1078         ),
1079         mesh,
1080         dimensionedScalar("cellId", dimless, 0.0)
1081     );
1082
1083     forAll(foamCellMap, cellI)
1084     {
1085         cellIdField[cellI] = foamCellMap[cellI];
1086     }
1087
1088     // Construct field with calculated bc to hold cell type.
1089     volScalarField cellTypeField
1090     (
1091         IOobject
1092         (
1093             "cellType",
1094             runTime.timeName(),
1095             mesh,
1096             IOobject::NO_READ,
1097             IOobject::AUTO_WRITE
1098         ),
1099         mesh,
1100         dimensionedScalar("cellType", dimless, 0.0)
1101     );
1102
1103     forAll(foamCellType, cellI)
1104     {
1105         cellTypeField[cellI] = foamCellType[cellI];
1106     }
1107
1108     Info<< "Writing cellIds as volScalarField to " << cellIdField.objectPath()
1109         << "..." << nl << endl;
1110     mesh.write();
1111
1112     Info<< "End\n" << endl;
1113
1114     return 0;
1115 }
1116
1117
1118 // ************************************************************************* //