applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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
  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     I-Deas unv format mesh conversion.
  26
  27     Uses either
  28     - DOF sets (757) or
  29     - face groups (2452(Cubit), 2467)
  30     to do patching.
  31     Works without but then puts all faces in defaultFaces patch.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "argList.H"
  36 #include "polyMesh.H"
  37 #include "Time.H"
  38 #include "IFstream.H"
  39 #include "cellModeller.H"
  40 #include "cellSet.H"
  41 #include "faceSet.H"
  42 #include "DynamicList.H"
  43
  44 #include <cassert>
  45 #include "MeshedSurfaces.H"
  46
  47 using namespace Foam;
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 namespace Foam
  52 {
  53     template<>
  54     inline unsigned Hash<face>::operator()(const face& t, unsigned seed) const
  55     {
  56         return Hasher(t.cdata(),t.size()*sizeof(label), seed);
  57     }
  58
  59     template<>
  60     inline unsigned Hash<face>::operator()(const face& t) const
  61     {
  62         return Hash<face>::operator()(t, 0);
  63     }
  64 }
  65 const string SEPARATOR("    -1");
  66
  67 bool isSeparator(const string& line)
  68 {
  69     return line.substr(0, 6) == SEPARATOR;
  70 }
  71
  72
  73 // Reads past -1 and reads element type
  74 label readTag(IFstream& is)
  75 {
  76     string tag;
  77     do
  78     {
  79         if (!is.good())
  80         {
  81             return -1;
  82         }
  83
  84         string line;
  85
  86         is.getLine(line);
  87
  88         if (line.size() < 6)
  89         {
  90             return -1;
  91         }
  92
  93         tag = line.substr(0, 6);
  94
  95     } while (tag == SEPARATOR);
  96
  97     return readLabel(IStringStream(tag)());
  98 }
  99
 100
 101 // Reads and prints header
 102 void readHeader(IFstream& is)
 103 {
 104     string line;
 105
 106     while (is.good())
 107     {
 108         is.getLine(line);
 109
 110         if (isSeparator(line))
 111         {
 112             break;
 113         }
 114         else
 115         {
 116             Info<< line << endl;
 117         }
 118     }
 119 }
 120
 121
 122 // Skip
 123 void skipSection(IFstream& is)
 124 {
 125     Info<< "Skipping section at line " << is.lineNumber() << '.' << endl;
 126
 127     string line;
 128
 129     while (is.good())
 130     {
 131         is.getLine(line);
 132
 133         if (isSeparator(line))
 134         {
 135             break;
 136         }
 137     }
 138 }
 139
 140
 141 scalar readUnvScalar(const string& unvString)
 142 {
 143     string s(unvString);
 144
 145     s.replaceAll("d", "E");
 146     s.replaceAll("D", "E");
 147
 148     return readScalar(IStringStream(s)());
 149 }
 150
 151
 152 // Reads unit section
 153 void readUnits
 154 (
 155     IFstream& is,
 156     scalar& lengthScale,
 157     scalar& forceScale,
 158     scalar& tempScale,
 159     scalar& tempOffset
 160 )
 161 {
 162     Info<< "Starting reading units at line " << is.lineNumber() << '.' << endl;
 163
 164     string line;
 165     is.getLine(line);
 166
 167     label l = readLabel(IStringStream(line.substr(0, 10))());
 168     Info<< "l:" << l << endl;
 169
 170     string units(line.substr(10, 20));
 171     Info<< "units:" << units << endl;
 172
 173     label unitType = readLabel(IStringStream(line.substr(30, 10))());
 174     Info<< "unitType:" << unitType << endl;
 175
 176     // Read lengthscales
 177     is.getLine(line);
 178
 179     lengthScale = readUnvScalar(line.substr(0, 25));
 180     forceScale = readUnvScalar(line.substr(25, 25));
 181     tempScale = readUnvScalar(line.substr(50, 25));
 182
 183     is.getLine(line);
 184     tempOffset = readUnvScalar(line.substr(0, 25));
 185
 186     Info<< "Unit factors:" << nl
 187         << "    Length scale       : " << lengthScale << nl
 188         << "    Force scale        : " << forceScale << nl
 189         << "    Temperature scale  : " << tempScale << nl
 190         << "    Temperature offset : " << tempOffset << nl
 191         << endl;
 192 }
 193
 194
 195 // Reads points section. Read region as well?
 196 void readPoints
 197 (
 198     IFstream& is,
 199     DynamicList<point>& points,     // coordinates
 200     DynamicList<label>& unvPointID  // unv index
 201 )
 202 {
 203     Info<< "Starting reading points at line " << is.lineNumber() << '.' << endl;
 204
 205     static bool hasWarned = false;
 206
 207     while (true)
 208     {
 209         string line;
 210         is.getLine(line);
 211
 212         label pointI = readLabel(IStringStream(line.substr(0, 10))());
 213
 214         if (pointI == -1)
 215         {
 216             break;
 217         }
 218         else if (pointI != points.size()+1 && !hasWarned)
 219         {
 220             hasWarned = true;
 221
 222             IOWarningIn
 223             (
 224                 "readPoints(IFstream&, label&, DynamicList<point>"
 225                 ", DynamicList<label>&)",
 226                 is
 227             )   << "Points not in order starting at point " << pointI
 228                 //<< " at line " << is.lineNumber()
 229                 //<< abort(FatalError);
 230                 << endl;
 231         }
 232
 233         point pt;
 234         is.getLine(line);
 235         pt[0] = readUnvScalar(line.substr(0, 25));
 236         pt[1] = readUnvScalar(line.substr(25, 25));
 237         pt[2] = readUnvScalar(line.substr(50, 25));
 238
 239         unvPointID.append(pointI);
 240         points.append(pt);
 241     }
 242
 243     points.shrink();
 244     unvPointID.shrink();
 245
 246     Info<< "Read " << points.size() << " points." << endl;
 247 }
 248
 249 void addAndExtend
 250 (
 251     DynamicList<label>& indizes,
 252     label cellI,
 253     label val
 254 )
 255 {
 256     if (indizes.size() < (cellI+1))
 257     {
 258         indizes.setSize(cellI+1,-1);
 259     }
 260     indizes[cellI] = val;
 261 }
 262
 263 // Reads cells section. Read region as well? Not handled yet but should just
 264 // be a matter of reading corresponding to boundaryFaces the correct property
 265 // and sorting it later on.
 266 void readCells
 267 (
 268     IFstream& is,
 269     DynamicList<cellShape>& cellVerts,
 270     DynamicList<label>& cellMaterial,
 271     DynamicList<label>& boundaryFaceIndices,
 272     DynamicList<face>& boundaryFaces,
 273     DynamicList<label>& cellCorrespondence,
 274     DynamicList<label>& unvPointID  // unv index
 275 )
 276 {
 277     Info<< "Starting reading cells at line " << is.lineNumber() << '.' << endl;
 278
 279     // Invert point numbering.
 280     label maxUnvPoint = 0;
 281     forAll(unvPointID, pointI)
 282     {
 283         maxUnvPoint = max(maxUnvPoint, unvPointID[pointI]);
 284     }
 285     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 286
 287
 288     const cellModel& hex = *(cellModeller::lookup("hex"));
 289     const cellModel& prism = *(cellModeller::lookup("prism"));
 290     const cellModel& tet = *(cellModeller::lookup("tet"));
 291
 292     labelHashSet skippedElements;
 293
 294     labelHashSet foundFeType;
 295
 296     while (true)
 297     {
 298         string line;
 299         is.getLine(line);
 300
 301         if (isSeparator(line))
 302         {
 303             break;
 304         }
 305
 306         label cellI, feID, physProp, matProp, colour, nNodes;
 307
 308         IStringStream lineStr(line);
 309         lineStr
 310             >> cellI >> feID >> physProp >> matProp >> colour >> nNodes;
 311
 312         if (foundFeType.insert(feID))
 313         {
 314             Info<< "First occurrence of element type " << feID
 315                 << " for cell " << cellI << " at line "
 316                 << is.lineNumber() << endl;
 317         }
 318
 319         if (feID == 11)
 320         {
 321             // Rod. Skip.
 322             is.getLine(line);
 323             is.getLine(line);
 324         }
 325         else if (feID == 171)
 326         {
 327             // Rod. Skip.
 328             is.getLine(line);
 329         }
 330         else if (feID == 41 || feID == 91)
 331         {
 332             // Triangle. Save - used for patching later on.
 333             is.getLine(line);
 334
 335             face cVerts(3);
 336             IStringStream lineStr(line);
 337             lineStr
 338                 >> cVerts[0] >> cVerts[1] >> cVerts[2];
 339             boundaryFaces.append(cVerts);
 340             boundaryFaceIndices.append(cellI);
 341         }
 342         else if (feID == 44 || feID == 94)
 343         {
 344             // Quad. Save - used for patching later on.
 345             is.getLine(line);
 346
 347             face cVerts(4);
 348             IStringStream lineStr(line);
 349             lineStr
 350                 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
 351             boundaryFaces.append(cVerts);
 352             boundaryFaceIndices.append(cellI);
 353         }
 354         else if (feID == 111)
 355         {
 356             // Tet.
 357             is.getLine(line);
 358
 359             labelList cVerts(4);
 360             IStringStream lineStr(line);
 361             lineStr
 362                 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
 363
 364             cellVerts.append(cellShape(tet, cVerts, true));
 365             cellMaterial.append(physProp);
 366             addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
 367
 368             if (cellVerts.last().size() != cVerts.size())
 369             {
 370                 Info<< "Line:" << is.lineNumber()
 371                     << " element:" << cellI
 372                     << " type:" << feID
 373                     << " collapsed from " << cVerts << nl
 374                     << " to:" << cellVerts.last()
 375                     << endl;
 376             }
 377         }
 378         else if (feID == 112)
 379         {
 380             // Wedge.
 381             is.getLine(line);
 382
 383             labelList cVerts(6);
 384             IStringStream lineStr(line);
 385             lineStr
 386                 >> cVerts[0] >> cVerts[1] >> cVerts[2]
 387                 >> cVerts[3] >> cVerts[4] >> cVerts[5];
 388
 389             cellVerts.append(cellShape(prism, cVerts, true));
 390             cellMaterial.append(physProp);
 391             addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
 392
 393             if (cellVerts.last().size() != cVerts.size())
 394             {
 395                 Info<< "Line:" << is.lineNumber()
 396                     << " element:" << cellI
 397                     << " type:" << feID
 398                     << " collapsed from " << cVerts << nl
 399                     << " to:" << cellVerts.last()
 400                     << endl;
 401             }
 402         }
 403         else if (feID == 115)
 404         {
 405             // Hex.
 406             is.getLine(line);
 407
 408             labelList cVerts(8);
 409             IStringStream lineStr(line);
 410             lineStr
 411                 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
 412                 >> cVerts[4] >> cVerts[5] >> cVerts[6] >> cVerts[7];
 413
 414             cellVerts.append(cellShape(hex, cVerts, true));
 415             cellMaterial.append(physProp);
 416             addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
 417
 418             if (cellVerts.last().size() != cVerts.size())
 419             {
 420                 Info<< "Line:" << is.lineNumber()
 421                     << " element:" << cellI
 422                     << " type:" << feID
 423                     << " collapsed from " << cVerts << nl
 424                     << " to:" << cellVerts.last()
 425                     << endl;
 426             }
 427         }
 428         else if (feID == 118)
 429         {
 430             // Parabolic Tet
 431             is.getLine(line);
 432
 433             labelList cVerts(4);
 434             label dummy;
 435             {
 436                 IStringStream lineStr(line);
 437                 lineStr
 438                     >> cVerts[0] >> dummy >> cVerts[1] >> dummy >> cVerts[2];
 439             }
 440             is.getLine(line);
 441             {
 442                 IStringStream lineStr(line);
 443                 lineStr >> dummy>> cVerts[3];
 444             }
 445
 446             cellVerts.append(cellShape(tet, cVerts, true));
 447             cellMaterial.append(physProp);
 448             addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
 449
 450             if (cellVerts.last().size() != cVerts.size())
 451             {
 452                 Info<< "Line:" << is.lineNumber()
 453                     << " element:" << cellI
 454                     << " type:" << feID
 455                     << " collapsed from " << cVerts << nl
 456                     << " to:" << cellVerts.last()
 457                     << endl;
 458             }
 459         }
 460         else
 461         {
 462             if (skippedElements.insert(feID))
 463             {
 464                 IOWarningIn("readCells(IFstream&, label&)", is)
 465                     << "Cell type " << feID << " not supported" << endl;
 466             }
 467             is.getLine(line);  // Do nothing
 468         }
 469     }
 470
 471     cellVerts.shrink();
 472     cellMaterial.shrink();
 473     boundaryFaces.shrink();
 474     boundaryFaceIndices.shrink();
 475     cellCorrespondence.shrink();
 476
 477     Info<< "Read " << cellVerts.size() << " cells"
 478         << " and " << boundaryFaces.size() << " boundary faces." << endl;
 479 }
 480
 481
 482 void readSets
 483 (
 484     IFstream& is,
 485     DynamicList<word>& patchNames,
 486     DynamicList<labelList>& patchFaceIndices
 487 )
 488 {
 489     Info<< "Starting reading patches at line " << is.lineNumber() << '.'
 490         << endl;
 491
 492     while (true)
 493     {
 494         string line;
 495         is.getLine(line);
 496
 497         if (isSeparator(line))
 498         {
 499             break;
 500         }
 501
 502         IStringStream lineStr(line);
 503         label group, constraintSet, restraintSet, loadSet, dofSet,
 504             tempSet, contactSet, nFaces;
 505         lineStr
 506             >> group >> constraintSet >> restraintSet >> loadSet
 507             >> dofSet >> tempSet >> contactSet >> nFaces;
 508
 509         is.getLine(line);
 510         word groupName = string::validate<word>(line);
 511
 512         Info<< "For group " << group
 513             << " named " << groupName
 514             << " trying to read " << nFaces << " patch face indices."
 515             << endl;
 516
 517         labelList groupIndices(nFaces);
 518         label groupType = -1;
 519         nFaces = 0;
 520
 521         while (nFaces < groupIndices.size())
 522         {
 523             is.getLine(line);
 524             IStringStream lineStr(line);
 525
 526             // Read one (for last face) or two entries from line.
 527             label nRead = 2;
 528             if (nFaces == groupIndices.size()-1)
 529             {
 530                 nRead = 1;
 531             }
 532
 533             for (label i = 0; i < nRead; i++)
 534             {
 535                 label tag, nodeLeaf, component;
 536
 537                 lineStr >> groupType >> tag >> nodeLeaf >> component;
 538
 539                 groupIndices[nFaces++] = tag;
 540             }
 541         }
 542
 543
 544         // Store
 545         if (groupType == 8)
 546         {
 547             patchNames.append(groupName);
 548             patchFaceIndices.append(groupIndices);
 549         }
 550         else
 551         {
 552             IOWarningIn("readSets(..)", is)
 553                 << "When reading patches expect entity type code 8"
 554                 << nl << "    Skipping group code " << groupType
 555                 << endl;
 556         }
 557     }
 558
 559     patchNames.shrink();
 560     patchFaceIndices.shrink();
 561 }
 562
 563
 564
 565 // Read dof set (757)
 566 void readDOFS
 567 (
 568     IFstream& is,
 569     DynamicList<word>& patchNames,
 570     DynamicList<labelList>& dofVertices
 571 )
 572 {
 573     Info<< "Starting reading contraints at line " << is.lineNumber() << '.'
 574         << endl;
 575
 576     string line;
 577     is.getLine(line);
 578     label group;
 579     {
 580         IStringStream lineStr(line);
 581         lineStr >> group;
 582     }
 583
 584     is.getLine(line);
 585     {
 586         IStringStream lineStr(line);
 587         patchNames.append(lineStr);
 588     }
 589
 590     Info<< "For DOF set " << group
 591         << " named " << patchNames.last()
 592         << " trying to read vertex indices."
 593         << endl;
 594
 595     DynamicList<label> vertices;
 596     while (true)
 597     {
 598         string line;
 599         is.getLine(line);
 600
 601         if (isSeparator(line))
 602         {
 603             break;
 604         }
 605
 606         IStringStream lineStr(line);
 607         label pointI;
 608         lineStr >> pointI;
 609
 610         vertices.append(pointI);
 611     }
 612
 613     Info<< "For DOF set " << group
 614         << " named " << patchNames.last()
 615         << " read " << vertices.size() << " vertex indices." << endl;
 616
 617     dofVertices.append(vertices.shrink());
 618
 619     patchNames.shrink();
 620     dofVertices.shrink();
 621 }
 622
 623
 624 // Returns -1 or group that all of the vertices of f are in,
 625 label findPatch(const List<labelHashSet>& dofGroups, const face& f)
 626 {
 627     forAll(dofGroups, patchI)
 628     {
 629         if (dofGroups[patchI].found(f[0]))
 630         {
 631             bool allInGroup = true;
 632
 633             // Check rest of face
 634             for (label fp = 1; fp < f.size(); fp++)
 635             {
 636                 if (!dofGroups[patchI].found(f[fp]))
 637                 {
 638                     allInGroup = false;
 639                     break;
 640                 }
 641             }
 642
 643             if (allInGroup)
 644             {
 645                 return patchI;
 646             }
 647         }
 648     }
 649     return -1;
 650 }
 651
 652
 653 // Main program:
 654
 655 int main(int argc, char *argv[])
 656 {
 657     argList::noParallel();
 658     argList::validArgs.append(".unv file");
 659     argList::addBoolOption
 660     (
 661         "dump",
 662         "dump boundary faces as boundaryFaces.obj (for debugging)"
 663     );
 664
 665 #   include "setRootCase.H"
 666 #   include "createTime.H"
 667
 668     const fileName ideasName = args[1];
 669     IFstream inFile(ideasName);
 670
 671     if (!inFile.good())
 672     {
 673         FatalErrorIn(args.executable())
 674             << "Cannot open file " << ideasName
 675             << exit(FatalError);
 676     }
 677
 678
 679     // Switch on additional debug info
 680     const bool verbose = false; //true;
 681
 682     // Unit scale factors
 683     scalar lengthScale = 1;
 684     scalar forceScale = 1;
 685     scalar tempScale = 1;
 686     scalar tempOffset = 0;
 687
 688     // Points
 689     DynamicList<point> points;
 690     // Original unv point label
 691     DynamicList<label> unvPointID;
 692
 693     // Cells
 694     DynamicList<cellShape> cellVerts;
 695     DynamicList<label> cellMat;
 696     DynamicList<label> cellCorrespondence;
 697
 698     // Boundary faces
 699     DynamicList<label> boundaryFaceIndices;
 700     DynamicList<face> boundaryFaces;
 701
 702     // Patch names and patchFace indices.
 703     DynamicList<word> patchNames;
 704     DynamicList<labelList> patchFaceIndices;
 705     DynamicList<labelList> dofVertIndices;
 706
 707
 708     while (inFile.good())
 709     {
 710         label tag = readTag(inFile);
 711
 712         if (tag == -1)
 713         {
 714             break;
 715         }
 716
 717         Info<< "Processing tag:" << tag << endl;
 718
 719         switch (tag)
 720         {
 721             case 151:
 722                 readHeader(inFile);
 723             break;
 724
 725             case 164:
 726                 readUnits
 727                 (
 728                     inFile,
 729                     lengthScale,
 730                     forceScale,
 731                     tempScale,
 732                     tempOffset
 733                 );
 734             break;
 735
 736             case 2411:
 737                 readPoints(inFile, points, unvPointID);
 738             break;
 739
 740             case 2412:
 741                 readCells
 742                 (
 743                     inFile,
 744                     cellVerts,
 745                     cellMat,
 746                     boundaryFaceIndices,
 747                     boundaryFaces,
 748                     cellCorrespondence,
 749                     unvPointID
 750                 );
 751             break;
 752
 753             case 2452:
 754             case 2467:
 755                 readSets
 756                 (
 757                     inFile,
 758                     patchNames,
 759                     patchFaceIndices
 760                 );
 761             break;
 762
 763             case 757:
 764                 readDOFS
 765                 (
 766                     inFile,
 767                     patchNames,
 768                     dofVertIndices
 769                 );
 770             break;
 771
 772             default:
 773                 Info<< "Skipping tag " << tag << " on line "
 774                     << inFile.lineNumber() << endl;
 775                 skipSection(inFile);
 776             break;
 777         }
 778         Info<< endl;
 779     }
 780
 781
 782     // Invert point numbering.
 783     label maxUnvPoint = 0;
 784     forAll(unvPointID, pointI)
 785     {
 786         maxUnvPoint = max(maxUnvPoint, unvPointID[pointI]);
 787     }
 788     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 789
 790
 791     // Renumber vertex numbers on cells
 792
 793     forAll(cellVerts, cellI)
 794     {
 795         labelList foamVerts
 796         (
 797             renumber
 798             (
 799                 unvToFoam,
 800                 static_cast<labelList&>(cellVerts[cellI])
 801             )
 802         );
 803
 804         if (findIndex(foamVerts, -1) != -1)
 805         {
 806             FatalErrorIn(args.executable())
 807                 << "Cell " << cellI
 808                 << " unv vertices " << cellVerts[cellI]
 809                 << " has some undefined vertices " << foamVerts
 810                 << abort(FatalError);
 811         }
 812
 813         // Bit nasty: replace vertex list.
 814         cellVerts[cellI].transfer(foamVerts);
 815     }
 816
 817     // Renumber vertex numbers on boundaryFaces
 818
 819     forAll(boundaryFaces, bFaceI)
 820     {
 821         labelList foamVerts(renumber(unvToFoam, boundaryFaces[bFaceI]));
 822
 823         if (findIndex(foamVerts, -1) != -1)
 824         {
 825             FatalErrorIn(args.executable())
 826                 << "Boundary face " << bFaceI
 827                 << " unv vertices " << boundaryFaces[bFaceI]
 828                 << " has some undefined vertices " << foamVerts
 829                 << abort(FatalError);
 830         }
 831
 832         // Bit nasty: replace vertex list.
 833         boundaryFaces[bFaceI].transfer(foamVerts);
 834     }
 835
 836
 837
 838     // Patchify = create patchFaceVerts for use in cellShape construction
 839     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 840
 841     // Works in one of two modes. Either has read boundaryFaces which
 842     // just need to be sorted according to patch. Or has read point constraint
 843     // sets (dofVertIndices).
 844
 845     List<faceList> patchFaceVerts;
 846
 847     labelList own(boundaryFaces.size(), -1);
 848     labelList nei(boundaryFaces.size(), -1);
 849     HashTable<label, label> faceToCell[2];
 850
 851     {
 852         HashTable<label, face, Hash<face> > faceToFaceID(boundaryFaces.size());
 853         forAll(boundaryFaces, faceI)
 854         {
 855             SortableList<label> sortedVerts(boundaryFaces[faceI]);
 856             faceToFaceID.insert(face(sortedVerts), faceI);
 857         }
 858
 859         forAll(cellVerts, cellI)
 860         {
 861             faceList faces = cellVerts[cellI].faces();
 862             forAll(faces, i)
 863             {
 864                 SortableList<label> sortedVerts(faces[i]);
 865                 HashTable<label, face, Hash<face> >::const_iterator fnd =
 866                     faceToFaceID.find(face(sortedVerts));
 867
 868                 if (fnd != faceToFaceID.end())
 869                 {
 870                     label faceI = fnd();
 871                     int stat = face::compare(faces[i], boundaryFaces[faceI]);
 872
 873                     if (stat == 1)
 874                     {
 875                         // Same orientation. Cell is owner.
 876                         own[faceI] = cellI;
 877                     }
 878                     else if (stat == -1)
 879                     {
 880                         // Opposite orientation. Cell is neighbour.
 881                         nei[faceI] = cellI;
 882                     }
 883                 }
 884             }
 885         }
 886
 887         label nReverse = 0;
 888         forAll(own, faceI)
 889         {
 890             if (own[faceI] == -1 && nei[faceI] != -1)
 891             {
 892                 // Boundary face with incorrect orientation
 893                 boundaryFaces[faceI] = boundaryFaces[faceI].reverseFace();
 894                 Swap(own[faceI], nei[faceI]);
 895                 nReverse++;
 896             }
 897         }
 898         if (nReverse > 0)
 899         {
 900             Info << "Found " << nReverse << " reversed boundary faces out of "
 901                 << boundaryFaces.size() << endl;
 902         }
 903
 904
 905         label cnt = 0;
 906         forAll(own, faceI)
 907         {
 908             if (own[faceI] != -1 && nei[faceI] != -1)
 909             {
 910                 cnt++;
 911             }
 912         }
 913
 914         if (cnt > 0)
 915         {
 916             Info << "Of " << boundaryFaces.size() << " so-called"
 917                 << " boundary faces " << cnt << " belong to two cells "
 918                 << "and are therefore internal" << endl;
 919         }
 920     }
 921
 922     HashTable<labelList,word> cellZones;
 923     HashTable<labelList,word> faceZones;
 924     List<bool> isAPatch(patchNames.size(),true);
 925
 926     if (dofVertIndices.size())
 927     {
 928         // Use the vertex constraints to patch. Is of course bit dodgy since
 929         // face goes on patch if all its vertices are on a constraint.
 930         // Note: very inefficient since goes through all faces (including
 931         // internal ones) twice. Should do a construct without patches
 932         // and then repatchify.
 933
 934         Info<< "Using " << dofVertIndices.size()
 935             << " DOF sets to detect boundary faces."<< endl;
 936
 937         // Renumber vertex numbers on contraints
 938         forAll(dofVertIndices, patchI)
 939         {
 940             inplaceRenumber(unvToFoam, dofVertIndices[patchI]);
 941         }
 942
 943
 944         // Build labelHashSet of points per dofGroup/patch
 945         List<labelHashSet> dofGroups(dofVertIndices.size());
 946
 947         forAll(dofVertIndices, patchI)
 948         {
 949             const labelList& foamVerts = dofVertIndices[patchI];
 950
 951             forAll(foamVerts, i)
 952             {
 953                 dofGroups[patchI].insert(foamVerts[i]);
 954             }
 955         }
 956
 957         List<DynamicList<face> > dynPatchFaces(dofVertIndices.size());
 958
 959         forAll(cellVerts, cellI)
 960         {
 961             const cellShape& shape = cellVerts[cellI];
 962
 963             const faceList shapeFaces(shape.faces());
 964
 965             forAll(shapeFaces, i)
 966             {
 967                 label patchI = findPatch(dofGroups, shapeFaces[i]);
 968
 969                 if (patchI != -1)
 970                 {
 971                     dynPatchFaces[patchI].append(shapeFaces[i]);
 972                 }
 973             }
 974         }
 975
 976         // Transfer
 977         patchFaceVerts.setSize(dynPatchFaces.size());
 978
 979         forAll(dynPatchFaces, patchI)
 980         {
 981             patchFaceVerts[patchI].transfer(dynPatchFaces[patchI]);
 982         }
 983     }
 984     else
 985     {
 986         // Use the boundary faces.
 987
 988         // Construct the patch faces by sorting the boundaryFaces according to
 989         // patch.
 990         patchFaceVerts.setSize(patchFaceIndices.size());
 991
 992         Info<< "Sorting boundary faces according to group (patch)" << endl;
 993
 994         // make sure that no face is used twice on the boundary
 995         // (possible for boundary-only faceZones)
 996         labelHashSet alreadyOnBoundary;
 997
 998         // Construct map from boundaryFaceIndices
 999         Map<label> boundaryFaceToIndex(boundaryFaceIndices.size());
1000
1001         forAll(boundaryFaceIndices, i)
1002         {
1003             boundaryFaceToIndex.insert(boundaryFaceIndices[i], i);
1004         }
1005
1006         forAll(patchFaceVerts, patchI)
1007         {
1008             Info << patchI << ": " << patchNames[patchI] << " is " << flush;
1009
1010             faceList& patchFaces = patchFaceVerts[patchI];
1011             const labelList& faceIndices = patchFaceIndices[patchI];
1012
1013             patchFaces.setSize(faceIndices.size());
1014
1015             bool duplicateFaces = false;
1016
1017             label cnt = 0;
1018             forAll(patchFaces, i)
1019             {
1020                 if (boundaryFaceToIndex.found(faceIndices[i]))
1021                 {
1022                     label bFaceI = boundaryFaceToIndex[faceIndices[i]];
1023
1024                     if (own[bFaceI] != -1 && nei[bFaceI] == -1)
1025                     {
1026                         patchFaces[cnt] = boundaryFaces[bFaceI];
1027                         cnt++;
1028                         if (alreadyOnBoundary.found(bFaceI))
1029                         {
1030                             duplicateFaces = true;
1031                         }
1032                     }
1033                 }
1034             }
1035
1036             if (cnt != patchFaces.size() || duplicateFaces)
1037             {
1038                 isAPatch[patchI] = false;
1039
1040                 if (verbose)
1041                 {
1042                     if (cnt != patchFaces.size())
1043                     {
1044                         WarningIn(args.executable())
1045                             << "For patch " << patchI << " there were "
1046                             << patchFaces.size()-cnt
1047                             << " faces not used because they seem"
1048                             << " to be internal. "
1049                             << "This seems to be a face or a cell-zone"
1050                             << endl;
1051                     }
1052                     else
1053                     {
1054                         WarningIn(args.executable())
1055                             << "Patch "
1056                             << patchI << " has faces that are already "
1057                             << " in use on other boundary-patches,"
1058                             << " Assuming faceZoneset." << endl;
1059                     }
1060                 }
1061
1062                 patchFaces.setSize(0); // Assume that this is no patch at all
1063
1064                 if (cellCorrespondence[faceIndices[0]] >= 0)
1065                 {
1066                     Info << "cellZone" << endl;
1067                     labelList theCells(faceIndices.size());
1068                     forAll(faceIndices, i)
1069                     {
1070                         if (cellCorrespondence[faceIndices[0]] < 0)
1071                         {
1072                             FatalErrorIn(args.executable())
1073                                 << "The face index " << faceIndices[i]
1074                                 << " was not found amongst the cells."
1075                                 << " This kills the theory that "
1076                                 << patchNames[patchI] << " is a cell zone"
1077                                 << endl
1078                                 << abort(FatalError);
1079                         }
1080                         theCells[i] = cellCorrespondence[faceIndices[i]];
1081                     }
1082                     cellZones.insert(patchNames[patchI], theCells);
1083                 }
1084                 else
1085                 {
1086                     Info << "faceZone" << endl;
1087                     labelList theFaces(faceIndices.size());
1088                     forAll(faceIndices, i)
1089                     {
1090                         theFaces[i] = boundaryFaceToIndex[faceIndices[i]];
1091                     }
1092                     faceZones.insert(patchNames[patchI],theFaces);
1093                 }
1094             }
1095             else
1096             {
1097                 Info << "patch" << endl;
1098
1099                 forAll(patchFaces, i)
1100                 {
1101                     label bFaceI = boundaryFaceToIndex[faceIndices[i]];
1102                     alreadyOnBoundary.insert(bFaceI);
1103                 }
1104             }
1105         }
1106     }
1107
1108     pointField polyPoints;
1109     polyPoints.transfer(points);
1110
1111     // Length scaling factor
1112     polyPoints /= lengthScale;
1113
1114
1115     // For debugging: dump boundary faces as OBJ surface mesh
1116     if (args.optionFound("dump"))
1117     {
1118         Info<< "Writing boundary faces to OBJ file boundaryFaces.obj"
1119             << nl << endl;
1120
1121         // Create globally numbered surface
1122         meshedSurface rawSurface
1123         (
1124             xferCopy(polyPoints),
1125             xferCopyTo< faceList >(boundaryFaces)
1126         );
1127
1128         // Write locally numbered surface
1129         meshedSurface
1130         (
1131             xferCopy(rawSurface.localPoints()),
1132             xferCopy(rawSurface.localFaces())
1133         ).write(runTime.path()/"boundaryFaces.obj");
1134     }
1135
1136
1137     Info<< "\nConstructing mesh with non-default patches of size:" << nl;
1138     DynamicList<word> usedPatchNames;
1139     DynamicList<faceList> usedPatchFaceVerts;
1140
1141     forAll(patchNames, patchI)
1142     {
1143         if (isAPatch[patchI])
1144         {
1145             Info<< "    " << patchNames[patchI] << '\t'
1146                 << patchFaceVerts[patchI].size() << nl;
1147             usedPatchNames.append(patchNames[patchI]);
1148             usedPatchFaceVerts.append(patchFaceVerts[patchI]);
1149         }
1150     }
1151     usedPatchNames.shrink();
1152     usedPatchFaceVerts.shrink();
1153
1154     Info<< endl;
1155
1156
1157
1158     // Construct mesh
1159     polyMesh mesh
1160     (
1161         IOobject
1162         (
1163             polyMesh::defaultRegion,
1164             runTime.constant(),
1165             runTime
1166         ),
1167         xferMove(polyPoints),
1168         cellVerts,
1169         usedPatchFaceVerts,             // boundaryFaces,
1170         usedPatchNames,                 // boundaryPatchNames,
1171         wordList(patchNames.size(), polyPatch::typeName), // boundaryPatchTypes,
1172         "defaultFaces",             // defaultFacesName
1173         polyPatch::typeName,        // defaultFacesType,
1174         wordList(0)                 // boundaryPatchPhysicalTypes
1175     );
1176
1177
1178     if (faceZones.size() > 0 || cellZones.size() > 0)
1179     {
1180         Info << "Adding cell and face zones" << endl;
1181
1182         List<pointZone*> pZones(0);
1183         List<faceZone*> fZones(faceZones.size());
1184         List<cellZone*> cZones(cellZones.size());
1185
1186         if (cellZones.size() > 0)
1187         {
1188             forAll(cellZones.toc(), cnt)
1189             {
1190                 word name = cellZones.toc()[cnt];
1191                 Info<< " Cell Zone " << name << " " << tab
1192                     << cellZones[name].size() << endl;
1193
1194                 cZones[cnt] = new cellZone
1195                 (
1196                     name,
1197                     cellZones[name],
1198                     cnt,
1199                     mesh.cellZones()
1200                 );
1201             }
1202         }
1203         if (faceZones.size() > 0)
1204         {
1205             const labelList& own = mesh.faceOwner();
1206             const labelList& nei = mesh.faceNeighbour();
1207             const pointField& centers = mesh.faceCentres();
1208             const pointField& points = mesh.points();
1209
1210             forAll(faceZones.toc(), cnt)
1211             {
1212                 word name = faceZones.toc()[cnt];
1213                 const labelList& oldIndizes = faceZones[name];
1214                 labelList indizes(oldIndizes.size());
1215
1216                 Info<< " Face Zone " << name << " " << tab
1217                     << oldIndizes.size() << endl;
1218
1219                 forAll(indizes, i)
1220                 {
1221                     const label old = oldIndizes[i];
1222                     label noveau = -1;
1223                     label c1 = -1, c2 = -1;
1224                     if (faceToCell[0].found(old))
1225                     {
1226                         c1 = faceToCell[0][old];
1227                     }
1228                     if (faceToCell[1].found(old))
1229                     {
1230                         c2 = faceToCell[1][old];
1231                     }
1232                     if (c1 < c2)
1233                     {
1234                         label tmp = c1;
1235                         c1 = c2;
1236                         c2 = tmp;
1237                     }
1238                     if (c2 == -1)
1239                     {
1240                         // Boundary face is part of the faceZone
1241                         forAll(own, j)
1242                         {
1243                             if (own[j] == c1)
1244                             {
1245                                 const face& f = boundaryFaces[old];
1246                                 if (mag(centers[j]- f.centre(points)) < SMALL)
1247                                 {
1248                                     noveau = j;
1249                                     break;
1250                                 }
1251                             }
1252                         }
1253                     }
1254                     else
1255                     {
1256                         forAll(nei, j)
1257                         {
1258                             if
1259                             (
1260                                 (c1 == own[j] && c2 == nei[j])
1261                              || (c2 == own[j] && c1 == nei[j])
1262                             )
1263                             {
1264                                 noveau = j;
1265                                 break;
1266                             }
1267                         }
1268                     }
1269                     assert(noveau > -1);
1270                     indizes[i] = noveau;
1271                 }
1272                 fZones[cnt] = new faceZone
1273                 (
1274                     faceZones.toc()[cnt],
1275                     indizes,
1276                     boolList(indizes.size(),false),
1277                     cnt,
1278                     mesh.faceZones()
1279                 );
1280             }
1281         }
1282         mesh.addZones(pZones, fZones, cZones);
1283
1284         Info << endl;
1285     }
1286
1287     mesh.write();
1288
1289     Info<< "End\n" << endl;
1290
1291     return 0;
1292 }
1293
1294
1295 // ************************************************************************* //