src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "processorPolyPatch.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "dictionary.H"
  29 #include "SubField.H"
  30 #include "demandDrivenData.H"
  31 #include "matchPoints.H"
  32 #include "OFstream.H"
  33 #include "polyMesh.H"
  34 #include "Time.H"
  35 #include "transformList.H"
  36 #include "PstreamBuffers.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(processorPolyPatch, 0);
  43     addToRunTimeSelectionTable(polyPatch, processorPolyPatch, dictionary);
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 Foam::processorPolyPatch::processorPolyPatch
  50 (
  51     const word& name,
  52     const label size,
  53     const label start,
  54     const label index,
  55     const polyBoundaryMesh& bm,
  56     const int myProcNo,
  57     const int neighbProcNo
  58 )
  59 :
  60     coupledPolyPatch(name, size, start, index, bm),
  61     myProcNo_(myProcNo),
  62     neighbProcNo_(neighbProcNo),
  63     neighbFaceCentres_(),
  64     neighbFaceAreas_(),
  65     neighbFaceCellCentres_()
  66 {}
  67
  68
  69 Foam::processorPolyPatch::processorPolyPatch
  70 (
  71     const word& name,
  72     const dictionary& dict,
  73     const label index,
  74     const polyBoundaryMesh& bm
  75 )
  76 :
  77     coupledPolyPatch(name, dict, index, bm),
  78     myProcNo_(readLabel(dict.lookup("myProcNo"))),
  79     neighbProcNo_(readLabel(dict.lookup("neighbProcNo"))),
  80     neighbFaceCentres_(),
  81     neighbFaceAreas_(),
  82     neighbFaceCellCentres_()
  83 {}
  84
  85
  86 Foam::processorPolyPatch::processorPolyPatch
  87 (
  88     const processorPolyPatch& pp,
  89     const polyBoundaryMesh& bm
  90 )
  91 :
  92     coupledPolyPatch(pp, bm),
  93     myProcNo_(pp.myProcNo_),
  94     neighbProcNo_(pp.neighbProcNo_),
  95     neighbFaceCentres_(),
  96     neighbFaceAreas_(),
  97     neighbFaceCellCentres_()
  98 {}
  99
 100
 101 Foam::processorPolyPatch::processorPolyPatch
 102 (
 103     const processorPolyPatch& pp,
 104     const polyBoundaryMesh& bm,
 105     const label index,
 106     const label newSize,
 107     const label newStart
 108 )
 109 :
 110     coupledPolyPatch(pp, bm, index, newSize, newStart),
 111     myProcNo_(pp.myProcNo_),
 112     neighbProcNo_(pp.neighbProcNo_),
 113     neighbFaceCentres_(),
 114     neighbFaceAreas_(),
 115     neighbFaceCellCentres_()
 116 {}
 117
 118
 119 Foam::processorPolyPatch::processorPolyPatch
 120 (
 121     const processorPolyPatch& pp,
 122     const polyBoundaryMesh& bm,
 123     const label index,
 124     const labelUList& mapAddressing,
 125     const label newStart
 126 )
 127 :
 128     coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
 129     myProcNo_(pp.myProcNo_),
 130     neighbProcNo_(pp.neighbProcNo_),
 131     neighbFaceCentres_(),
 132     neighbFaceAreas_(),
 133     neighbFaceCellCentres_()
 134 {}
 135
 136
 137 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 138
 139 Foam::processorPolyPatch::~processorPolyPatch()
 140 {
 141     neighbPointsPtr_.clear();
 142     neighbEdgesPtr_.clear();
 143 }
 144
 145
 146 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 147
 148 void Foam::processorPolyPatch::initGeometry(PstreamBuffers& pBufs)
 149 {
 150     if (Pstream::parRun())
 151     {
 152         UOPstream toNeighbProc(neighbProcNo(), pBufs);
 153
 154         toNeighbProc
 155             << faceCentres()
 156             << faceAreas()
 157             << faceCellCentres();
 158     }
 159 }
 160
 161
 162 void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 163 {
 164     if (Pstream::parRun())
 165     {
 166         {
 167             UIPstream fromNeighbProc(neighbProcNo(), pBufs);
 168
 169             fromNeighbProc
 170                 >> neighbFaceCentres_
 171                 >> neighbFaceAreas_
 172                 >> neighbFaceCellCentres_;
 173         }
 174
 175
 176         // My normals
 177         vectorField faceNormals(size());
 178
 179         // Neighbour normals
 180         vectorField nbrFaceNormals(neighbFaceAreas_.size());
 181
 182         // Calculate normals from areas and check
 183         forAll(faceNormals, facei)
 184         {
 185             scalar magSf = mag(faceAreas()[facei]);
 186             scalar nbrMagSf = mag(neighbFaceAreas_[facei]);
 187             scalar avSf = (magSf + nbrMagSf)/2.0;
 188
 189             if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
 190             {
 191                 // Undetermined normal. Use dummy normal to force separation
 192                 // check. (note use of sqrt(VSMALL) since that is how mag
 193                 // scales)
 194                 faceNormals[facei] = point(1, 0, 0);
 195                 nbrFaceNormals[facei] = faceNormals[facei];
 196             }
 197             else if (mag(magSf - nbrMagSf)/avSf > matchTolerance())
 198             {
 199                 fileName nm
 200                 (
 201                     boundaryMesh().mesh().time().path()
 202                    /name()+"_faces.obj"
 203                 );
 204                 Pout<< "processorPolyPatch::order : Writing my " << size()
 205                     << " faces to OBJ file " << nm << endl;
 206                 writeOBJ(nm, *this, points());
 207
 208                 FatalErrorIn
 209                 (
 210                     "processorPolyPatch::calcGeometry()"
 211                 )   << "face " << facei << " area does not match neighbour by "
 212                     << 100*mag(magSf - nbrMagSf)/avSf
 213                     << "% -- possible face ordering problem." << endl
 214                     << "patch:" << name()
 215                     << " my area:" << magSf
 216                     << " neighbour area:" << nbrMagSf
 217                     << " matching tolerance:" << matchTolerance()
 218                     << endl
 219                     << "Mesh face:" << start()+facei
 220                     << " vertices:"
 221                     << UIndirectList<point>(points(), operator[](facei))()
 222                     << endl
 223                     << "If you are certain your matching is correct"
 224                     << " you can increase the 'matchTolerance' setting"
 225                     << " in the patch dictionary in the boundary file."
 226                     << endl
 227                     << "Rerun with processor debug flag set for"
 228                     << " more information." << exit(FatalError);
 229             }
 230             else
 231             {
 232                 faceNormals[facei] = faceAreas()[facei]/magSf;
 233                 nbrFaceNormals[facei] = neighbFaceAreas_[facei]/nbrMagSf;
 234             }
 235         }
 236
 237         calcTransformTensors
 238         (
 239             faceCentres(),
 240             neighbFaceCentres_,
 241             faceNormals,
 242             nbrFaceNormals,
 243             calcFaceTol(matchTolerance(), *this, points(), faceCentres()),
 244             matchTolerance()
 245         );
 246     }
 247 }
 248
 249
 250 void Foam::processorPolyPatch::initMovePoints
 251 (
 252     PstreamBuffers& pBufs,
 253     const pointField& p
 254 )
 255 {
 256     polyPatch::movePoints(pBufs, p);
 257     processorPolyPatch::initGeometry(pBufs);
 258 }
 259
 260
 261 void Foam::processorPolyPatch::movePoints
 262 (
 263     PstreamBuffers& pBufs,
 264     const pointField&
 265 )
 266 {
 267     processorPolyPatch::calcGeometry(pBufs);
 268 }
 269
 270
 271 void Foam::processorPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 272 {
 273     polyPatch::initUpdateMesh(pBufs);
 274
 275     if (Pstream::parRun())
 276     {
 277         // Express all points as patch face and index in face.
 278         labelList pointFace(nPoints());
 279         labelList pointIndex(nPoints());
 280
 281         for (label patchPointI = 0; patchPointI < nPoints(); patchPointI++)
 282         {
 283             label faceI = pointFaces()[patchPointI][0];
 284
 285             pointFace[patchPointI] = faceI;
 286
 287             const face& f = localFaces()[faceI];
 288
 289             pointIndex[patchPointI] = findIndex(f, patchPointI);
 290         }
 291
 292         // Express all edges as patch face and index in face.
 293         labelList edgeFace(nEdges());
 294         labelList edgeIndex(nEdges());
 295
 296         for (label patchEdgeI = 0; patchEdgeI < nEdges(); patchEdgeI++)
 297         {
 298             label faceI = edgeFaces()[patchEdgeI][0];
 299
 300             edgeFace[patchEdgeI] = faceI;
 301
 302             const labelList& fEdges = faceEdges()[faceI];
 303
 304             edgeIndex[patchEdgeI] = findIndex(fEdges, patchEdgeI);
 305         }
 306
 307         UOPstream toNeighbProc(neighbProcNo(), pBufs);
 308
 309         toNeighbProc
 310             << pointFace
 311             << pointIndex
 312             << edgeFace
 313             << edgeIndex;
 314     }
 315 }
 316
 317
 318 void Foam::processorPolyPatch::updateMesh(PstreamBuffers& pBufs)
 319 {
 320     // For completeness
 321     polyPatch::updateMesh(pBufs);
 322
 323     neighbPointsPtr_.clear();
 324     neighbEdgesPtr_.clear();
 325
 326     if (Pstream::parRun())
 327     {
 328         labelList nbrPointFace;
 329         labelList nbrPointIndex;
 330         labelList nbrEdgeFace;
 331         labelList nbrEdgeIndex;
 332
 333         {
 334             // Note cannot predict exact size since opposite nPoints might
 335             // be different from one over here.
 336             UIPstream fromNeighbProc(neighbProcNo(), pBufs);
 337
 338             fromNeighbProc
 339                 >> nbrPointFace
 340                 >> nbrPointIndex
 341                 >> nbrEdgeFace
 342                 >> nbrEdgeIndex;
 343         }
 344
 345         // Convert neighbour faces and indices into face back into
 346         // my edges and points.
 347
 348         // Convert points.
 349         // ~~~~~~~~~~~~~~~
 350
 351         neighbPointsPtr_.reset(new labelList(nPoints(), -1));
 352         labelList& neighbPoints = neighbPointsPtr_();
 353
 354         forAll(nbrPointFace, nbrPointI)
 355         {
 356             // Find face and index in face on this side.
 357             const face& f = localFaces()[nbrPointFace[nbrPointI]];
 358             label index = (f.size() - nbrPointIndex[nbrPointI]) % f.size();
 359             label patchPointI = f[index];
 360
 361             if (neighbPoints[patchPointI] == -1)
 362             {
 363                 // First reference of point
 364                 neighbPoints[patchPointI] = nbrPointI;
 365             }
 366             else if (neighbPoints[patchPointI] >= 0)
 367             {
 368                 // Point already visited. Mark as duplicate.
 369                 neighbPoints[patchPointI] = -2;
 370             }
 371         }
 372
 373         // Reset all duplicate entries to -1.
 374         forAll(neighbPoints, patchPointI)
 375         {
 376             if (neighbPoints[patchPointI] == -2)
 377             {
 378                 neighbPoints[patchPointI] = -1;
 379             }
 380         }
 381
 382         // Convert edges.
 383         // ~~~~~~~~~~~~~~
 384
 385         neighbEdgesPtr_.reset(new labelList(nEdges(), -1));
 386         labelList& neighbEdges = neighbEdgesPtr_();
 387
 388         forAll(nbrEdgeFace, nbrEdgeI)
 389         {
 390             // Find face and index in face on this side.
 391             const labelList& f = faceEdges()[nbrEdgeFace[nbrEdgeI]];
 392             label index = (f.size() - nbrEdgeIndex[nbrEdgeI] - 1) % f.size();
 393             label patchEdgeI = f[index];
 394
 395             if (neighbEdges[patchEdgeI] == -1)
 396             {
 397                 // First reference of edge
 398                 neighbEdges[patchEdgeI] = nbrEdgeI;
 399             }
 400             else if (neighbEdges[patchEdgeI] >= 0)
 401             {
 402                 // Edge already visited. Mark as duplicate.
 403                 neighbEdges[patchEdgeI] = -2;
 404             }
 405         }
 406
 407         // Reset all duplicate entries to -1.
 408         forAll(neighbEdges, patchEdgeI)
 409         {
 410             if (neighbEdges[patchEdgeI] == -2)
 411             {
 412                 neighbEdges[patchEdgeI] = -1;
 413             }
 414         }
 415
 416         // Remove any addressing used for shared points/edges calculation
 417         // since mostly not needed.
 418         primitivePatch::clearOut();
 419     }
 420 }
 421
 422
 423 const Foam::labelList& Foam::processorPolyPatch::neighbPoints() const
 424 {
 425     if (!neighbPointsPtr_.valid())
 426     {
 427         FatalErrorIn("processorPolyPatch::neighbPoints() const")
 428             << "No extended addressing calculated for patch " << name()
 429             << abort(FatalError);
 430     }
 431     return neighbPointsPtr_();
 432 }
 433
 434
 435 const Foam::labelList& Foam::processorPolyPatch::neighbEdges() const
 436 {
 437     if (!neighbEdgesPtr_.valid())
 438     {
 439         FatalErrorIn("processorPolyPatch::neighbEdges() const")
 440             << "No extended addressing calculated for patch " << name()
 441             << abort(FatalError);
 442     }
 443     return neighbEdgesPtr_();
 444 }
 445
 446
 447 void Foam::processorPolyPatch::initOrder
 448 (
 449     PstreamBuffers& pBufs,
 450     const primitivePatch& pp
 451 ) const
 452 {
 453     if (!Pstream::parRun())
 454     {
 455         return;
 456     }
 457
 458     if (debug)
 459     {
 460         fileName nm
 461         (
 462             boundaryMesh().mesh().time().path()
 463            /name()+"_faces.obj"
 464         );
 465         Pout<< "processorPolyPatch::order : Writing my " << pp.size()
 466             << " faces to OBJ file " << nm << endl;
 467         writeOBJ(nm, pp, pp.points());
 468
 469         // Calculate my face centres
 470         const pointField& fc = pp.faceCentres();
 471
 472         OFstream localStr
 473         (
 474             boundaryMesh().mesh().time().path()
 475            /name() + "_localFaceCentres.obj"
 476         );
 477         Pout<< "processorPolyPatch::order : "
 478             << "Dumping " << fc.size()
 479             << " local faceCentres to " << localStr.name() << endl;
 480
 481         forAll(fc, faceI)
 482         {
 483             writeOBJ(localStr, fc[faceI]);
 484         }
 485     }
 486
 487     if (owner())
 488     {
 489         pointField anchors(getAnchorPoints(pp, pp.points()));
 490
 491         // Now send all info over to the neighbour
 492         UOPstream toNeighbour(neighbProcNo(), pBufs);
 493         toNeighbour << pp.faceCentres() << anchors;
 494     }
 495 }
 496
 497
 498 // Return new ordering. Ordering is -faceMap: for every face index
 499 // the new face -rotation:for every new face the clockwise shift
 500 // of the original face. Return false if nothing changes (faceMap
 501 // is identity, rotation is 0)
 502 bool Foam::processorPolyPatch::order
 503 (
 504     PstreamBuffers& pBufs,
 505     const primitivePatch& pp,
 506     labelList& faceMap,
 507     labelList& rotation
 508 ) const
 509 {
 510     // Note: we only get the faces that originate from internal faces.
 511
 512     if (!Pstream::parRun())
 513     {
 514         return false;
 515     }
 516
 517     faceMap.setSize(pp.size());
 518     faceMap = -1;
 519
 520     rotation.setSize(pp.size());
 521     rotation = 0;
 522
 523     if (owner())
 524     {
 525         // Do nothing (i.e. identical mapping, zero rotation).
 526         // See comment at top.
 527         forAll(faceMap, patchFaceI)
 528         {
 529             faceMap[patchFaceI] = patchFaceI;
 530         }
 531
 532         return false;
 533     }
 534     else
 535     {
 536         vectorField masterCtrs;
 537         vectorField masterAnchors;
 538
 539         // Receive data from neighbour
 540         {
 541             UIPstream fromNeighbour(neighbProcNo(), pBufs);
 542             fromNeighbour >> masterCtrs >> masterAnchors;
 543         }
 544
 545         // Calculate typical distance from face centre
 546         scalarField tols
 547         (
 548             calcFaceTol(matchTolerance(), pp, pp.points(), pp.faceCentres())
 549         );
 550
 551         if (debug || masterCtrs.size() != pp.size())
 552         {
 553             {
 554                 OFstream nbrStr
 555                 (
 556                     boundaryMesh().mesh().time().path()
 557                    /name() + "_nbrFaceCentres.obj"
 558                 );
 559                 Pout<< "processorPolyPatch::order : "
 560                     << "Dumping neighbour faceCentres to " << nbrStr.name()
 561                     << endl;
 562                 forAll(masterCtrs, faceI)
 563                 {
 564                     writeOBJ(nbrStr, masterCtrs[faceI]);
 565                 }
 566             }
 567
 568             if (masterCtrs.size() != pp.size())
 569             {
 570                 FatalErrorIn
 571                 (
 572                     "processorPolyPatch::order(const primitivePatch&"
 573                     ", labelList&, labelList&) const"
 574                 )   << "in patch:" << name() << " : "
 575                     << "Local size of patch is " << pp.size() << " (faces)."
 576                     << endl
 577                     << "Received from neighbour " << masterCtrs.size()
 578                     << " faceCentres!"
 579                     << abort(FatalError);
 580             }
 581         }
 582
 583         // Geometric match of face centre vectors
 584         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 585
 586         // 1. Try existing ordering and transformation
 587         bool matchedAll = matchPoints
 588         (
 589             pp.faceCentres(),
 590             masterCtrs,
 591             tols,
 592             true,
 593             faceMap
 594         );
 595
 596         if (!matchedAll || debug)
 597         {
 598             // Dump faces
 599             fileName str
 600             (
 601                 boundaryMesh().mesh().time().path()
 602                /name()/name()+"_faces.obj"
 603             );
 604             Pout<< "processorPolyPatch::order :"
 605                 << " Writing faces to OBJ file " << str.name() << endl;
 606             writeOBJ(str, pp, pp.points());
 607
 608             OFstream ccStr
 609             (
 610                 boundaryMesh().mesh().time().path()
 611                /name() + "_faceCentresConnections.obj"
 612             );
 613
 614             Pout<< "processorPolyPatch::order :"
 615                 << " Dumping newly found match as lines between"
 616                 << " corresponding face centres to OBJ file " << ccStr.name()
 617                 << endl;
 618
 619             label vertI = 0;
 620
 621             forAll(pp.faceCentres(), faceI)
 622             {
 623                 label masterFaceI = faceMap[faceI];
 624
 625                 if (masterFaceI != -1)
 626                 {
 627                     const point& c0 = masterCtrs[masterFaceI];
 628                     const point& c1 = pp.faceCentres()[faceI];
 629                     writeOBJ(ccStr, c0, c1, vertI);
 630                 }
 631             }
 632         }
 633
 634         if (!matchedAll)
 635         {
 636             SeriousErrorIn
 637             (
 638                 "processorPolyPatch::order(const primitivePatch&"
 639                 ", labelList&, labelList&) const"
 640             )   << "in patch:" << name() << " : "
 641                 << "Cannot match vectors to faces on both sides of patch"
 642                 << endl
 643                 << "    masterCtrs[0]:" << masterCtrs[0] << endl
 644                 << "    ctrs[0]:" << pp.faceCentres()[0] << endl
 645                 << "    Please check your topology changes or maybe you have"
 646                 << " multiple separated (from cyclics) processor patches"
 647                 << endl
 648                 << "    Continuing with incorrect face ordering from now on!"
 649                 << endl;
 650
 651             return false;
 652         }
 653
 654         // Set rotation.
 655         forAll(faceMap, oldFaceI)
 656         {
 657             // The face f will be at newFaceI (after morphing) and we want its
 658             // anchorPoint (= f[0]) to align with the anchorpoint for the
 659             // corresponding face on the other side.
 660
 661             label newFaceI = faceMap[oldFaceI];
 662
 663             const point& wantedAnchor = masterAnchors[newFaceI];
 664
 665             rotation[newFaceI] = getRotation
 666             (
 667                 pp.points(),
 668                 pp[oldFaceI],
 669                 wantedAnchor,
 670                 tols[oldFaceI]
 671             );
 672
 673             if (rotation[newFaceI] == -1)
 674             {
 675                 SeriousErrorIn
 676                 (
 677                     "processorPolyPatch::order(const primitivePatch&"
 678                     ", labelList&, labelList&) const"
 679                 )   << "in patch " << name()
 680                     << " : "
 681                     << "Cannot find point on face " << pp[oldFaceI]
 682                     << " with vertices "
 683                     << UIndirectList<point>(pp.points(), pp[oldFaceI])()
 684                     << " that matches point " << wantedAnchor
 685                     << " when matching the halves of processor patch " << name()
 686                     << "Continuing with incorrect face ordering from now on!"
 687                     << endl;
 688
 689                 return false;
 690             }
 691         }
 692
 693         forAll(faceMap, faceI)
 694         {
 695             if (faceMap[faceI] != faceI || rotation[faceI] != 0)
 696             {
 697                 return true;
 698             }
 699         }
 700
 701         return false;
 702     }
 703 }
 704
 705
 706 void Foam::processorPolyPatch::write(Ostream& os) const
 707 {
 708     coupledPolyPatch::write(os);
 709     os.writeKeyword("myProcNo") << myProcNo_
 710         << token::END_STATEMENT << nl;
 711     os.writeKeyword("neighbProcNo") << neighbProcNo_
 712         << token::END_STATEMENT << nl;
 713 }
 714
 715
 716 // ************************************************************************* //