applications/solvers/surfaceTracking/freeSurface/freeSurfacePointDisplacement.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 \*---------------------------------------------------------------------------*/
  25
  26
  27 #include "freeSurface.H"
  28 #include "primitivePatchInterpolation.H"
  29 #include "emptyFaPatch.H"
  30 #include "wedgeFaPatch.H"
  31 #include "PstreamCombineReduceOps.H"
  32 #include "coordinateSystem.H"
  33 #include "scalarMatrices.H"
  34
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43
  44 tmp<vectorField> freeSurface::pointDisplacement(const scalarField& deltaH)
  45 {
  46     const pointField& points = aMesh().patch().localPoints();
  47     const labelListList& pointFaces = aMesh().patch().pointFaces();
  48
  49     controlPoints() += facesDisplacementDir()*deltaH;
  50
  51     tmp<vectorField> tdisplacement
  52     (
  53         new vectorField
  54         (
  55             points.size(),
  56             vector::zero
  57         )
  58     );
  59
  60     vectorField& displacement = tdisplacement();
  61
  62
  63     // Calculate displacement of internal points
  64     const vectorField& pointNormals = aMesh().pointAreaNormals();
  65     const edgeList& edges = aMesh().patch().edges();
  66     labelList internalPoints = aMesh().internalPoints();
  67
  68     forAll (internalPoints, pointI)
  69     {
  70         label curPoint = internalPoints[pointI];
  71
  72         const labelList& curPointFaces = pointFaces[curPoint];
  73
  74         vectorField lsPoints(curPointFaces.size(), vector::zero);
  75
  76         for (label i=0; i<curPointFaces.size(); i++)
  77         {
  78             label curFace = curPointFaces[i];
  79
  80             lsPoints[i] = controlPoints()[curFace];
  81         }
  82
  83         vectorField pointAndNormal =
  84             lsPlanePointAndNormal
  85             (
  86                 lsPoints,
  87                 points[curPoint],
  88                 pointNormals[curPoint]
  89             );
  90
  91         vector& P = pointAndNormal[0];
  92         vector& N = pointAndNormal[1];
  93
  94         displacement[curPoint] =
  95             pointsDisplacementDir()[curPoint]
  96            *((P - points[curPoint])&N)
  97            /(pointsDisplacementDir()[curPoint]&N);
  98     }
  99
 100
 101     // Mirror control points
 102     FieldField<Field, vector> patchMirrorPoints(aMesh().boundary().size());
 103
 104     forAll(patchMirrorPoints, patchI)
 105     {
 106         patchMirrorPoints.set
 107         (
 108             patchI,
 109             new vectorField
 110             (
 111                 aMesh().boundary()[patchI].faPatch::size(),
 112                 vector::zero
 113             )
 114         );
 115
 116         vectorField N =
 117             aMesh().boundary()[patchI].ngbPolyPatchFaceNormals();
 118
 119         const labelList peFaces =
 120             labelList::subList
 121             (
 122                 aMesh().edgeOwner(),
 123                 aMesh().boundary()[patchI].faPatch::size(),
 124                 aMesh().boundary()[patchI].start()
 125             );
 126
 127         const labelList& pEdges = aMesh().boundary()[patchI];
 128
 129         vectorField peCentres(pEdges.size(), vector::zero);
 130         forAll(peCentres, edgeI)
 131         {
 132             peCentres[edgeI] =
 133                 edges[pEdges[edgeI]].centre(points);
 134         }
 135
 136         vectorField delta =
 137             vectorField(controlPoints(), peFaces)
 138           - peCentres;
 139
 140         patchMirrorPoints[patchI] =
 141             peCentres + ((I - 2*N*N)&delta);
 142     }
 143
 144
 145     // Calculate displacement of boundary points
 146     labelList boundaryPoints = aMesh().boundaryPoints();
 147
 148     const labelListList& edgeFaces = aMesh().patch().edgeFaces();
 149     const labelListList& pointEdges = aMesh().patch().pointEdges();
 150
 151     forAll (boundaryPoints, pointI)
 152     {
 153         label curPoint = boundaryPoints[pointI];
 154
 155         if (motionPointsMask()[curPoint] == 1)
 156         {
 157             // Calculating mirror points
 158             const labelList& curPointEdges = pointEdges[curPoint];
 159
 160             vectorField mirrorPoints(2, vector::zero);
 161
 162             label counter = -1;
 163
 164             forAll (curPointEdges, edgeI)
 165             {
 166                 label curEdge = curPointEdges[edgeI];
 167
 168                 if(edgeFaces[curEdge].size() == 1)
 169                 {
 170                     label patchID = -1;
 171                     label edgeID = -1;
 172                     forAll(aMesh().boundary(), patchI)
 173                     {
 174                         const labelList& pEdges =
 175                             aMesh().boundary()[patchI];
 176                         label index = findIndex(pEdges, curEdge);
 177                         if (index != -1)
 178                         {
 179                             patchID = patchI;
 180                             edgeID = index;
 181                             break;
 182                         }
 183                     }
 184
 185                     mirrorPoints[++counter] =
 186                         patchMirrorPoints[patchID][edgeID];
 187                 }
 188             }
 189
 190             // Calculating LS plane fit
 191             const labelList& curPointFaces = pointFaces[curPoint];
 192
 193             vectorField lsPoints
 194             (
 195                 curPointFaces.size() + mirrorPoints.size(),
 196                 vector::zero
 197             );
 198
 199             counter = -1;
 200
 201             for (label i=0; i<curPointFaces.size(); i++)
 202             {
 203                 label curFace = curPointFaces[i];
 204
 205                 lsPoints[++counter] = controlPoints()[curFace];
 206             }
 207
 208             for (label i=0; i<mirrorPoints.size(); i++)
 209             {
 210                 lsPoints[++counter] = mirrorPoints[i];
 211             }
 212
 213             vectorField pointAndNormal =
 214                 lsPlanePointAndNormal
 215                 (
 216                     lsPoints,
 217                     points[curPoint],
 218                     pointNormals[curPoint]
 219                 );
 220
 221             vector& P = pointAndNormal[0];
 222             vector& N = pointAndNormal[1];
 223
 224             displacement[curPoint] =
 225                 pointsDisplacementDir()[curPoint]
 226                *((P - points[curPoint])&N)
 227                /(pointsDisplacementDir()[curPoint]&N);
 228         }
 229     }
 230
 231
 232     // Calculate displacement of axis point
 233     forAll (aMesh().boundary(), patchI)
 234     {
 235         if
 236         (
 237             aMesh().boundary()[patchI].type()
 238          == wedgeFaPatch::typeName
 239         )
 240         {
 241             const wedgeFaPatch& wedgePatch =
 242                 refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 243
 244             if(wedgePatch.axisPoint() > -1)
 245             {
 246                 label axisPoint = wedgePatch.axisPoint();
 247
 248                 displacement[axisPoint] =
 249                     pointsDisplacementDir()[axisPoint]
 250                    *(
 251                         pointsDisplacementDir()[axisPoint]
 252                        &(
 253                             controlPoints()[pointFaces[axisPoint][0]]
 254                           - points[axisPoint]
 255                         )
 256                     );
 257             }
 258         }
 259     }
 260
 261
 262     // Calculate displacement of processor patch points
 263     forAll (aMesh().boundary(), patchI)
 264     {
 265         if
 266         (
 267             aMesh().boundary()[patchI].type()
 268          == processorFaPatch::typeName
 269         )
 270         {
 271             const processorFaPatch& procPatch =
 272                 refCast<const processorFaPatch>(aMesh().boundary()[patchI]);
 273
 274             const labelList& patchPointLabels =
 275                 procPatch.pointLabels();
 276
 277             FieldField<Field, vector> lsPoints(patchPointLabels.size());
 278             forAll(lsPoints, pointI)
 279             {
 280                 lsPoints.set(pointI, new vectorField(0, vector::zero));
 281             }
 282
 283             const labelList& nonGlobalPatchPoints =
 284                 procPatch.nonGlobalPatchPoints();
 285
 286             forAll(nonGlobalPatchPoints, pointI)
 287             {
 288                 label curPatchPoint =
 289                     nonGlobalPatchPoints[pointI];
 290
 291                 label curPoint =
 292                     patchPointLabels[curPatchPoint];
 293
 294                 const labelList& curPointFaces = pointFaces[curPoint];
 295
 296                 lsPoints[curPatchPoint].setSize(curPointFaces.size());
 297
 298                 forAll(curPointFaces, faceI)
 299                 {
 300                     label curFace = curPointFaces[faceI];
 301
 302                     lsPoints[curPatchPoint][faceI] = controlPoints()[curFace];
 303                 }
 304
 305 #               include "boundaryProcessorFaPatchPoints.H"
 306             }
 307
 308             scalar lsPointsSize = 0;
 309             forAll(lsPoints, pointI)
 310             {
 311                 lsPointsSize +=
 312                     2*lsPoints[pointI].size()*sizeof(vector);
 313             }
 314
 315             // Parallel data exchange
 316             {
 317                 OPstream toNeighbProc
 318                 (
 319                     Pstream::blocking,
 320                     procPatch.neighbProcNo(),
 321                     lsPointsSize
 322                 );
 323
 324                 toNeighbProc << lsPoints;
 325             }
 326
 327             FieldField<Field, vector> ngbLsPoints(patchPointLabels.size());
 328
 329             {
 330                 IPstream fromNeighbProc
 331                 (
 332                     Pstream::blocking,
 333                     procPatch.neighbProcNo(),
 334                     lsPointsSize
 335                 );
 336
 337                 fromNeighbProc >> ngbLsPoints;
 338             }
 339
 340             forAll(nonGlobalPatchPoints, pointI)
 341             {
 342                 label curPatchPoint =
 343                     nonGlobalPatchPoints[pointI];
 344
 345                 label curPoint =
 346                     patchPointLabels[curPatchPoint];
 347
 348                 label curNgbPoint = procPatch.neighbPoints()[curPatchPoint];
 349
 350                 vectorField allLsPoints
 351                 (
 352                     lsPoints[curPatchPoint].size()
 353                   + ngbLsPoints[curNgbPoint].size(),
 354                     vector::zero
 355                 );
 356
 357                 label counter = -1;
 358                 forAll(lsPoints[curPatchPoint], pointI)
 359                 {
 360                     allLsPoints[++counter] = lsPoints[curPatchPoint][pointI];
 361                 }
 362                 forAll(ngbLsPoints[curNgbPoint], pointI)
 363                 {
 364                     allLsPoints[++counter] = ngbLsPoints[curNgbPoint][pointI];
 365                 }
 366
 367                 vectorField pointAndNormal =
 368                     lsPlanePointAndNormal
 369                     (
 370                         allLsPoints,
 371                         points[curPoint],
 372                         pointNormals[curPoint]
 373                     );
 374
 375                 vector& P = pointAndNormal[0];
 376                 vector& N = pointAndNormal[1];
 377
 378                 if (motionPointsMask()[curPoint] != 0)
 379                 {
 380                     displacement[curPoint] =
 381                         pointsDisplacementDir()[curPoint]
 382                        *((P - points[curPoint])&N)
 383                        /(pointsDisplacementDir()[curPoint]&N);
 384                 }
 385             }
 386         }
 387     }
 388
 389
 390     // Calculate displacement of global processor patch points
 391     if (aMesh().globalData().nGlobalPoints() > 0)
 392     {
 393         const labelList& spLabels =
 394             aMesh().globalData().sharedPointLabels();
 395
 396         const labelList& addr = aMesh().globalData().sharedPointAddr();
 397
 398         for (label k=0; k<aMesh().globalData().nGlobalPoints(); k++)
 399         {
 400             List<List<vector> > procLsPoints(Pstream::nProcs());
 401
 402             label curSharedPointIndex = findIndex(addr, k);
 403
 404             if (curSharedPointIndex != -1)
 405             {
 406                 label curPoint = spLabels[curSharedPointIndex];
 407
 408                 const labelList& curPointFaces = pointFaces[curPoint];
 409
 410                 procLsPoints[Pstream::myProcNo()] =
 411                     List<vector>(curPointFaces.size());
 412
 413                 forAll (curPointFaces, faceI)
 414                 {
 415                     label curFace = curPointFaces[faceI];
 416
 417                     procLsPoints[Pstream::myProcNo()][faceI] =
 418                         controlPoints()[curFace];
 419                 }
 420             }
 421
 422             Pstream::gatherList(procLsPoints);
 423             Pstream::scatterList(procLsPoints);
 424
 425             if (curSharedPointIndex != -1)
 426             {
 427                 label curPoint = spLabels[curSharedPointIndex];
 428
 429                 label nAllPoints = 0;
 430                 forAll(procLsPoints, procI)
 431                 {
 432                     nAllPoints += procLsPoints[procI].size();
 433                 }
 434
 435                 vectorField allPoints(nAllPoints, vector::zero);
 436
 437                 label counter = 0;
 438                 forAll(procLsPoints, procI)
 439                 {
 440                     forAll(procLsPoints[procI], pointI)
 441                     {
 442                         allPoints[counter++] =
 443                             procLsPoints[procI][pointI];
 444                     }
 445                 }
 446
 447                 vectorField pointAndNormal =
 448                     lsPlanePointAndNormal
 449                     (
 450                         allPoints,
 451                         points[curPoint],
 452                         pointNormals[curPoint]
 453                     );
 454
 455                 const vector& P = pointAndNormal[0];
 456                 const vector& N = pointAndNormal[1];
 457
 458                 displacement[curPoint] =
 459                     pointsDisplacementDir()[curPoint]
 460                    *((P - points[curPoint])&N)
 461                    /(pointsDisplacementDir()[curPoint]&N);
 462             }
 463         }
 464     }
 465
 466     return tdisplacement;
 467 }
 468
 469
 470 tmp<vectorField> freeSurface::lsPlanePointAndNormal
 471 (
 472     const vectorField& points,
 473     const vector& origin,
 474     const vector& axis
 475 ) const
 476 {
 477     // LS in local CS
 478     vector dir = (points[0] - origin);
 479     dir -= axis*(axis&dir);
 480     dir /= mag(dir);
 481     coordinateSystem cs("cs", origin, axis, dir);
 482
 483     vectorField localPoints = cs.localPosition(points);
 484     scalarField W = 1.0/(mag(points - origin) + SMALL);
 485
 486     scalarRectangularMatrix M
 487     (
 488         points.size(),
 489         3,
 490         0.0
 491     );
 492
 493     for (label i=0; i<localPoints.size(); i++)
 494     {
 495         M[i][0] = localPoints[i].x();
 496         M[i][1] = localPoints[i].y();
 497         M[i][2] = 1.0;
 498     }
 499
 500     scalarSquareMatrix MtM(3, 0.0);
 501     for (label i = 0; i < MtM.n(); i++)
 502     {
 503         for (label j = 0; j < MtM.m(); j++)
 504         {
 505             for (label k = 0; k < M.n(); k++)
 506             {
 507                 MtM[i][j] += M[k][i]*M[k][j]*W[k];
 508             }
 509         }
 510     }
 511
 512     scalarField MtR(3, 0);
 513     for (label i = 0; i < MtR.size(); i++)
 514     {
 515         for (label j = 0; j < M.n(); j++)
 516         {
 517             MtR[i] += M[j][i]*localPoints[j].z()*W[j];
 518         }
 519     }
 520
 521     scalarSquareMatrix::LUsolve(MtM, MtR);
 522
 523     vector n0 = vector(-MtR[0], -MtR[1], 1);
 524     n0 = cs.globalVector(n0);
 525     n0 /= mag(n0);
 526
 527     vector p0 = vector(0, 0, MtR[2]);
 528     p0 = cs.globalPosition(p0);
 529
 530     tmp<vectorField> pointAndNormal
 531     (
 532         new vectorField(2, vector::zero)
 533     );
 534
 535     pointAndNormal()[0] = p0;
 536     pointAndNormal()[1] = n0;
 537
 538     return pointAndNormal;
 539 }
 540
 541
 542 // tmp<vectorField> freeSurface::pointDisplacementSM()
 543 // {
 544 //     const pointField& points = aMesh().patch().localPoints();
 545 //     const labelListList& pointFaces = aMesh().patch().pointFaces();
 546
 547
 548 //     tmp<vectorField> tdisplacement
 549 //     (
 550 //         new vectorField
 551 //         (
 552 //             points.size(),
 553 //             vector::zero
 554 //         )
 555 //     );
 556
 557 //     vectorField& displacement = tdisplacement();
 558
 559
 560 //     forAll (pointFaces, pointI)
 561 //     {
 562 //         scalar weightsSum = 0.0;
 563
 564 //         const labelList& curPointFaces = pointFaces[pointI];
 565
 566 //         forAll (curPointFaces, faceI)
 567 //         {
 568 //             label curFace = curPointFaces[faceI];
 569
 570 //             scalar weight = 1.0/mag
 571 //             (
 572 //                 points[pointI]
 573 //               - controlPoints()[curFace]
 574 //             );
 575
 576 //             displacement[pointI] += weight*controlPoints()[curFace];
 577
 578 //             weightsSum += weight;
 579 //         }
 580
 581 //         displacement[pointI] /= weightsSum;
 582
 583 //         displacement[pointI] -= points[pointI];
 584 //     }
 585
 586
 587 //     displacement = motionPointsMask()*
 588 //         (pointsDisplacementDir()&displacement)*
 589 //         pointsDisplacementDir();
 590
 591
 592 //     return tdisplacement;
 593 // }
 594
 595
 596 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 597
 598 } // End namespace Foam
 599
 600 // ************************************************************************* //