src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "interactionLists.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 Foam::scalar Foam::interactionLists::transTol = 1e-12;
  32
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 void Foam::interactionLists::buildCellReferralLists()
  37 {
  38     Info<< nl << "Determining molecule referring schedule" << endl;
  39
  40     const referredCellList& refIntL(ril());
  41
  42     DynamicList<label> referralProcs;
  43
  44     // Run through all referredCells to build list of interacting processors
  45
  46     forAll(refIntL, rIL)
  47     {
  48         const referredCell& rC(refIntL[rIL]);
  49
  50         if (findIndex(referralProcs, rC.sourceProc()) == -1)
  51         {
  52             referralProcs.append(rC.sourceProc());
  53         }
  54     }
  55
  56     referralProcs.shrink();
  57
  58 //     Pout << "referralProcs: " << nl << referralProcs << endl;
  59
  60     List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
  61
  62     List<DynamicList<DynamicList<label> > >
  63         cellReceivingReferralLists(referralProcs.size());
  64
  65     // Run through all referredCells again building up send and receive info
  66
  67     forAll(refIntL, rIL)
  68     {
  69         const referredCell& rC(refIntL[rIL]);
  70
  71         label rPI = findIndex(referralProcs, rC.sourceProc());
  72
  73         DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
  74
  75         DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
  76
  77         label existingSource = findIndex(sRL, rC.sourceCell());
  78
  79         // Check to see if this source cell has already been allocated to
  80         // come to this processor.  If not, add the source cell to the sending
  81         // list and add the current referred cell to the receiving list.
  82
  83         // It shouldn't be possible for the sending and receiving lists to be
  84         // different lengths, because their append operations happen at the
  85         // same time.
  86
  87         if (existingSource == -1)
  88         {
  89             sRL.append(rC.sourceCell());
  90
  91             rRL.append
  92             (
  93                 DynamicList<label> (labelList(1,rIL))
  94             );
  95         }
  96         else
  97         {
  98             rRL[existingSource].append(rIL);
  99
 100             rRL[existingSource].shrink();
 101         }
 102     }
 103
 104     forAll(referralProcs, rPI)
 105     {
 106         DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
 107
 108         DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
 109
 110         sRL.shrink();
 111
 112         rRL.shrink();
 113     }
 114
 115     // It is assumed that cell exchange is reciprocal, if proc A has cells to
 116     // send to proc B, then proc B must have some to send to proc A.
 117
 118     cellReceivingReferralLists_.setSize(referralProcs.size());
 119
 120     cellSendingReferralLists_.setSize(referralProcs.size());
 121
 122     forAll(referralProcs, rPI)
 123     {
 124         DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
 125
 126         labelListList translLL(rRL.size());
 127
 128         forAll(rRL, rRLI)
 129         {
 130             translLL[rRLI] = rRL[rRLI];
 131         }
 132
 133         cellReceivingReferralLists_[rPI] = receivingReferralList
 134         (
 135             referralProcs[rPI],
 136             translLL
 137         );
 138     }
 139
 140     // Send sendingReferralLists to each interacting processor.
 141
 142     forAll(referralProcs, rPI)
 143     {
 144
 145         DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
 146
 147         if (referralProcs[rPI] != Pstream::myProcNo())
 148         {
 149             if (Pstream::parRun())
 150             {
 151                 OPstream toInteractingProc
 152                 (
 153                     Pstream::blocking,
 154                     referralProcs[rPI]
 155                 );
 156
 157                 toInteractingProc << sendingReferralList
 158                 (
 159                     Pstream::myProcNo(),
 160                     sRL
 161                 );
 162             }
 163         }
 164     }
 165
 166     // Receive sendingReferralLists from each interacting processor.
 167
 168     forAll(referralProcs, rPI)
 169     {
 170         if (referralProcs[rPI] != Pstream::myProcNo())
 171         {
 172             if (Pstream::parRun())
 173             {
 174                 IPstream fromInteractingProc
 175                 (
 176                     Pstream::blocking,
 177                     referralProcs[rPI]
 178                 );
 179
 180                 fromInteractingProc >> cellSendingReferralLists_[rPI];
 181             }
 182         }
 183         else
 184         {
 185             cellSendingReferralLists_[rPI] = sendingReferralList
 186             (
 187                 Pstream::myProcNo(),
 188                 cellSendingReferralLists[rPI]
 189             );
 190         }
 191     }
 192 }
 193
 194
 195 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 196
 197
 198 Foam::interactionLists::interactionLists
 199 (
 200     const polyMesh& mesh,
 201     scalar rCutMaxSqr,
 202     bool pointPointListBuild
 203 )
 204 :
 205     mesh_(mesh),
 206     rCutMaxSqr_(rCutMaxSqr),
 207     dil_(*this, pointPointListBuild),
 208     ril_(*this, pointPointListBuild),
 209     cellSendingReferralLists_(),
 210     cellReceivingReferralLists_()
 211 {
 212     buildCellReferralLists();
 213 }
 214
 215
 216 Foam::interactionLists::interactionLists(const polyMesh& mesh)
 217 :
 218     mesh_(mesh),
 219     dil_(*this),
 220     ril_(*this)
 221 {}
 222
 223
 224 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 225
 226 Foam::interactionLists::~interactionLists()
 227 {}
 228
 229
 230 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 231
 232 bool Foam::interactionLists::testPointPointDistance
 233 (
 234     const label ptI,
 235     const label ptJ
 236 ) const
 237 {
 238     return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_);
 239 }
 240
 241
 242 bool Foam::interactionLists::testEdgeEdgeDistance
 243 (
 244     const edge& eI,
 245     const edge& eJ
 246 ) const
 247 {
 248     const vector& eJs(mesh_.points()[eJ.start()]);
 249     const vector& eJe(mesh_.points()[eJ.end()]);
 250
 251     return testEdgeEdgeDistance(eI, eJs, eJe);
 252 }
 253
 254
 255 bool Foam::interactionLists::testPointFaceDistance
 256 (
 257     const label p,
 258     const label faceNo
 259 ) const
 260 {
 261     const vector& pointPosition(mesh_.points()[p]);
 262
 263     return testPointFaceDistance(pointPosition, faceNo);
 264 }
 265
 266
 267 bool Foam::interactionLists::testPointFaceDistance
 268 (
 269     const label p,
 270     const referredCell& refCell
 271 ) const
 272 {
 273     const vector& pointPosition(mesh_.points()[p]);
 274
 275     forAll (refCell.faces(), rCF)
 276     {
 277         if
 278         (
 279             testPointFaceDistance
 280             (
 281                 pointPosition,
 282                 refCell.faces()[rCF],
 283                 refCell.vertexPositions(),
 284                 refCell.faceCentres()[rCF],
 285                 refCell.faceAreas()[rCF]
 286             )
 287         )
 288         {
 289             return true;
 290         }
 291     }
 292
 293     return false;
 294 }
 295
 296
 297 bool Foam::interactionLists::testPointFaceDistance
 298 (
 299     const vectorList& pointsToTest,
 300     const label faceNo
 301 ) const
 302 {
 303     forAll(pointsToTest, pTT)
 304     {
 305         const vector& p(pointsToTest[pTT]);
 306
 307         // if any point in the list is in range of the face
 308         // then the rest do not need to be tested and
 309         // true can be returned
 310
 311         if (testPointFaceDistance(p, faceNo))
 312         {
 313             return true;
 314         }
 315     }
 316
 317     return false;
 318 }
 319
 320
 321 bool Foam::interactionLists::testPointFaceDistance
 322 (
 323     const vector& p,
 324     const label faceNo
 325 ) const
 326 {
 327     const face& faceToTest(mesh_.faces()[faceNo]);
 328
 329     const vector& faceC(mesh_.faceCentres()[faceNo]);
 330
 331     const vector& faceA(mesh_.faceAreas()[faceNo]);
 332
 333     const vectorList& points(mesh_.points());
 334
 335     return testPointFaceDistance
 336     (
 337         p,
 338         faceToTest,
 339         points,
 340         faceC,
 341         faceA
 342     );
 343 }
 344
 345
 346 bool Foam::interactionLists::testPointFaceDistance
 347 (
 348     const vector& p,
 349     const labelList& faceToTest,
 350     const vectorList& points,
 351     const vector& faceC,
 352     const vector& faceA
 353 ) const
 354 {
 355     vector faceN(faceA/mag(faceA));
 356
 357     scalar perpDist((p - faceC) & faceN);
 358
 359     if (magSqr(perpDist) > rCutMaxSqr_)
 360     {
 361         return false;
 362     }
 363
 364     vector pointOnPlane = (p - faceN * perpDist);
 365
 366     if (magSqr(faceC - pointOnPlane) < rCutMaxSqr_*1e-8)
 367     {
 368         // If pointOnPlane is very close to the face centre
 369         // then defining the local axes will be inaccurate
 370         // and it is very likely that pointOnPlane will be
 371         // inside the face, so return true if the points
 372         // are in range to be safe
 373
 374         return (magSqr(pointOnPlane - p) <= rCutMaxSqr_);
 375     }
 376
 377     vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
 378
 379     vector yAxis =
 380         ((faceC - pointOnPlane) ^ faceN)
 381        /mag((faceC - pointOnPlane) ^ faceN);
 382
 383     List<vector2D> local2DVertices(faceToTest.size());
 384
 385     forAll(faceToTest, fTT)
 386     {
 387         const vector& V(points[faceToTest[fTT]]);
 388
 389         if (magSqr(V-p) <= rCutMaxSqr_)
 390         {
 391             return true;
 392         }
 393
 394         local2DVertices[fTT] = vector2D
 395         (
 396             ((V - pointOnPlane) & xAxis),
 397             ((V - pointOnPlane) & yAxis)
 398         );
 399     }
 400
 401     scalar localFaceCx((faceC - pointOnPlane) & xAxis);
 402
 403     scalar la_valid = -1;
 404
 405     forAll(local2DVertices, fV)
 406     {
 407         const vector2D& va(local2DVertices[fV]);
 408
 409         const vector2D& vb
 410         (
 411             local2DVertices[(fV + 1) % local2DVertices.size()]
 412         );
 413
 414         if (mag(vb.y()-va.y()) > SMALL)
 415         {
 416             scalar la =
 417                 (
 418                     va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
 419                 )
 420                /localFaceCx;
 421
 422             scalar lv = -va.y()/(vb.y() - va.y());
 423
 424
 425             if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
 426             {
 427                 la_valid = la;
 428
 429                 break;
 430             }
 431         }
 432     }
 433
 434     if (la_valid < 0)
 435     {
 436         // perpendicular point inside face, nearest point is pointOnPlane;
 437         return (magSqr(pointOnPlane-p) <= rCutMaxSqr_);
 438     }
 439     else
 440     {
 441         // perpendicular point outside face, nearest point is
 442         // on edge that generated la_valid;
 443         return
 444         (
 445             magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
 446          <= rCutMaxSqr_
 447         );
 448     }
 449
 450     // if the algorithm hasn't returned anything by now then something has
 451     // gone wrong.
 452
 453     FatalErrorIn("interactionLists.C") << nl
 454         << "point " << p << " to face " << faceToTest
 455         << " comparison did not find a nearest point"
 456         << " to be inside or outside face."
 457         << abort(FatalError);
 458
 459     return false;
 460 }
 461
 462
 463 bool Foam::interactionLists::testEdgeEdgeDistance
 464 (
 465     const edge& eI,
 466     const vector& eJs,
 467     const vector& eJe
 468 ) const
 469 {
 470     vector a(eI.vec(mesh_.points()));
 471     vector b(eJe - eJs);
 472
 473     const vector& eIs(mesh_.points()[eI.start()]);
 474
 475     vector c(eJs - eIs);
 476
 477     vector crossab = a ^ b;
 478     scalar magCrossSqr = magSqr(crossab);
 479
 480     if (magCrossSqr > VSMALL)
 481     {
 482         // If the edges are parallel then a point-face
 483         // search will pick them up
 484
 485         scalar s = ((c ^ b) & crossab)/magCrossSqr;
 486         scalar t = ((c ^ a) & crossab)/magCrossSqr;
 487
 488         // Check for end points outside of range 0..1
 489         // If the closest point is outside this range
 490         // a point-face search will have found it.
 491
 492         return
 493         (
 494             s >= 0
 495          && s <= 1
 496          && t >= 0
 497          && t <= 1
 498          && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr_
 499         );
 500     }
 501
 502     return false;
 503 }
 504
 505
 506 const Foam::labelList Foam::interactionLists::realCellsInRangeOfSegment
 507 (
 508     const labelList& segmentFaces,
 509     const labelList& segmentEdges,
 510     const labelList& segmentPoints
 511 ) const
 512 {
 513     DynamicList<label> realCellsFoundInRange;
 514
 515     forAll(segmentFaces, sF)
 516     {
 517         const label f = segmentFaces[sF];
 518
 519         forAll (mesh_.points(), p)
 520         {
 521             if (testPointFaceDistance(p, f))
 522             {
 523                 const labelList& pCells(mesh_.pointCells()[p]);
 524
 525                 forAll(pCells, pC)
 526                 {
 527                     const label cellI(pCells[pC]);
 528
 529                     if (findIndex(realCellsFoundInRange, cellI) == -1)
 530                     {
 531                         realCellsFoundInRange.append(cellI);
 532                     }
 533                 }
 534             }
 535         }
 536     }
 537
 538     forAll(segmentPoints, sP)
 539     {
 540         const label p = segmentPoints[sP];
 541
 542         forAll(mesh_.faces(), f)
 543         {
 544             if (testPointFaceDistance(p, f))
 545             {
 546                 const label cellO(mesh_.faceOwner()[f]);
 547
 548                 if (findIndex(realCellsFoundInRange, cellO) == -1)
 549                 {
 550                     realCellsFoundInRange.append(cellO);
 551                 }
 552
 553                 if (mesh_.isInternalFace(f))
 554                 {
 555                     // boundary faces will not have neighbour information
 556
 557                     const label cellN(mesh_.faceNeighbour()[f]);
 558
 559                     if (findIndex(realCellsFoundInRange, cellN) == -1)
 560                     {
 561                         realCellsFoundInRange.append(cellN);
 562                     }
 563                 }
 564             }
 565         }
 566     }
 567
 568     forAll(segmentEdges, sE)
 569     {
 570         const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
 571
 572         forAll (mesh_.edges(), edgeIIndex)
 573         {
 574             const edge& eI(mesh_.edges()[edgeIIndex]);
 575
 576             if (testEdgeEdgeDistance(eI, eJ))
 577             {
 578                 const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
 579
 580                 forAll(eICells, eIC)
 581                 {
 582                     const label cellI(eICells[eIC]);
 583
 584                     if (findIndex(realCellsFoundInRange, cellI) == -1)
 585                     {
 586                         realCellsFoundInRange.append(cellI);
 587                     }
 588                 }
 589             }
 590         }
 591     }
 592
 593     return realCellsFoundInRange.shrink();
 594 }
 595
 596
 597 const Foam::labelList Foam::interactionLists::referredCellsInRangeOfSegment
 598 (
 599     const List<referredCell>& referredInteractionList,
 600     const labelList& segmentFaces,
 601     const labelList& segmentEdges,
 602     const labelList& segmentPoints
 603 ) const
 604 {
 605     DynamicList<label> referredCellsFoundInRange;
 606
 607     forAll(segmentFaces, sF)
 608     {
 609         const label f = segmentFaces[sF];
 610
 611         forAll(referredInteractionList, rIL)
 612         {
 613             const vectorList& refCellPoints
 614                 = referredInteractionList[rIL].vertexPositions();
 615
 616             if (testPointFaceDistance(refCellPoints, f))
 617             {
 618                 if (findIndex(referredCellsFoundInRange, rIL) == -1)
 619                 {
 620                     referredCellsFoundInRange.append(rIL);
 621                 }
 622             }
 623         }
 624     }
 625
 626     forAll(segmentPoints, sP)
 627     {
 628         const label p = segmentPoints[sP];
 629
 630         forAll(referredInteractionList, rIL)
 631         {
 632             const referredCell& refCell(referredInteractionList[rIL]);
 633
 634             if (testPointFaceDistance(p, refCell))
 635             {
 636                 if (findIndex(referredCellsFoundInRange, rIL) == -1)
 637                 {
 638                     referredCellsFoundInRange.append(rIL);
 639                 }
 640             }
 641         }
 642     }
 643
 644     forAll(segmentEdges, sE)
 645     {
 646         const edge& eI(mesh_.edges()[segmentEdges[sE]]);
 647
 648         forAll(referredInteractionList, rIL)
 649         {
 650             const vectorList& refCellPoints
 651                 = referredInteractionList[rIL].vertexPositions();
 652
 653             const edgeList& refCellEdges
 654                 = referredInteractionList[rIL].edges();
 655
 656             forAll(refCellEdges, rCE)
 657             {
 658                 const edge& eJ(refCellEdges[rCE]);
 659
 660                 if
 661                 (
 662                     testEdgeEdgeDistance
 663                     (
 664                         eI,
 665                         refCellPoints[eJ.start()],
 666                         refCellPoints[eJ.end()]
 667                     )
 668                 )
 669                 {
 670                     if(findIndex(referredCellsFoundInRange, rIL) == -1)
 671                     {
 672                         referredCellsFoundInRange.append(rIL);
 673                     }
 674                 }
 675             }
 676         }
 677     }
 678
 679     return referredCellsFoundInRange.shrink();
 680 }
 681
 682
 683 // ************************************************************************* //