src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.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 "isoSurface.H"
  27 #include "polyMesh.H"
  28 #include "syncTools.H"
  29 #include "surfaceFields.H"
  30 #include "OFstream.H"
  31 #include "meshTools.H"
  32
  33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  34
  35 template<class Type>
  36 Foam::tmp<Foam::SlicedGeometricField
  37 <
  38     Type,
  39     Foam::fvPatchField,
  40     Foam::slicedFvPatchField,
  41     Foam::volMesh
  42 > >
  43 Foam::isoSurface::adaptPatchFields
  44 (
  45     const GeometricField<Type, fvPatchField, volMesh>& fld
  46 ) const
  47 {
  48     typedef SlicedGeometricField
  49     <
  50         Type,
  51         fvPatchField,
  52         slicedFvPatchField,
  53         volMesh
  54     > FieldType;
  55
  56     tmp<FieldType> tsliceFld
  57     (
  58         new FieldType
  59         (
  60             IOobject
  61             (
  62                 fld.name(),
  63                 fld.instance(),
  64                 fld.db(),
  65                 IOobject::NO_READ,
  66                 IOobject::NO_WRITE,
  67                 false
  68             ),
  69             fld,        // internal field
  70             true        // preserveCouples
  71         )
  72     );
  73     FieldType& sliceFld = tsliceFld();
  74
  75     const fvMesh& mesh = fld.mesh();
  76
  77     const polyBoundaryMesh& patches = mesh.boundaryMesh();
  78
  79     forAll(patches, patchI)
  80     {
  81         const polyPatch& pp = patches[patchI];
  82
  83         if
  84         (
  85             isA<emptyPolyPatch>(pp)
  86          && pp.size() != sliceFld.boundaryField()[patchI].size()
  87         )
  88         {
  89             // Clear old value. Cannot resize it since is a slice.
  90             sliceFld.boundaryField().set(patchI, NULL);
  91
  92             // Set new value we can change
  93             sliceFld.boundaryField().set
  94             (
  95                 patchI,
  96                 new calculatedFvPatchField<Type>
  97                 (
  98                     mesh.boundary()[patchI],
  99                     sliceFld
 100                 )
 101             );
 102
 103             // Note: cannot use patchInternalField since uses emptyFvPatch::size
 104             // Do our own internalField instead.
 105             const labelUList& faceCells =
 106                 mesh.boundary()[patchI].patch().faceCells();
 107
 108             Field<Type>& pfld = sliceFld.boundaryField()[patchI];
 109             pfld.setSize(faceCells.size());
 110             forAll(faceCells, i)
 111             {
 112                 pfld[i] = sliceFld[faceCells[i]];
 113             }
 114         }
 115         else if (isA<cyclicPolyPatch>(pp))
 116         {
 117             // Already has interpolate as value
 118         }
 119         else if (isA<processorPolyPatch>(pp))
 120         {
 121             fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
 122             (
 123                 fld.boundaryField()[patchI]
 124             );
 125
 126             const scalarField& w = mesh.weights().boundaryField()[patchI];
 127
 128             tmp<Field<Type> > f =
 129                 w*pfld.patchInternalField()
 130               + (1.0-w)*pfld.patchNeighbourField();
 131
 132             PackedBoolList isCollocated
 133             (
 134                 collocatedFaces(refCast<const processorPolyPatch>(pp))
 135             );
 136
 137             forAll(isCollocated, i)
 138             {
 139                 if (!isCollocated[i])
 140                 {
 141                     pfld[i] = f()[i];
 142                 }
 143             }
 144         }
 145     }
 146     return tsliceFld;
 147 }
 148
 149
 150 template<class Type>
 151 Type Foam::isoSurface::generatePoint
 152 (
 153     const scalar s0,
 154     const Type& p0,
 155     const bool hasSnap0,
 156     const Type& snapP0,
 157
 158     const scalar s1,
 159     const Type& p1,
 160     const bool hasSnap1,
 161     const Type& snapP1
 162 ) const
 163 {
 164     scalar d = s1-s0;
 165
 166     if (mag(d) > VSMALL)
 167     {
 168         scalar s = (iso_-s0)/d;
 169
 170         if (hasSnap1 && s >= 0.5 && s <= 1)
 171         {
 172             return snapP1;
 173         }
 174         else if (hasSnap0 && s >= 0.0 && s <= 0.5)
 175         {
 176             return snapP0;
 177         }
 178         else
 179         {
 180             return s*p1 + (1.0-s)*p0;
 181         }
 182     }
 183     else
 184     {
 185         scalar s = 0.4999;
 186
 187         return s*p1 + (1.0-s)*p0;
 188     }
 189 }
 190
 191
 192 template<class Type>
 193 void Foam::isoSurface::generateTriPoints
 194 (
 195     const scalar s0,
 196     const Type& p0,
 197     const bool hasSnap0,
 198     const Type& snapP0,
 199
 200     const scalar s1,
 201     const Type& p1,
 202     const bool hasSnap1,
 203     const Type& snapP1,
 204
 205     const scalar s2,
 206     const Type& p2,
 207     const bool hasSnap2,
 208     const Type& snapP2,
 209
 210     const scalar s3,
 211     const Type& p3,
 212     const bool hasSnap3,
 213     const Type& snapP3,
 214
 215     DynamicList<Type>& points
 216 ) const
 217 {
 218     int triIndex = 0;
 219     if (s0 < iso_)
 220     {
 221         triIndex |= 1;
 222     }
 223     if (s1 < iso_)
 224     {
 225         triIndex |= 2;
 226     }
 227     if (s2 < iso_)
 228     {
 229         triIndex |= 4;
 230     }
 231     if (s3 < iso_)
 232     {
 233         triIndex |= 8;
 234     }
 235
 236     /* Form the vertices of the triangles for each case */
 237     switch (triIndex)
 238     {
 239         case 0x00:
 240         case 0x0F:
 241         break;
 242
 243         case 0x0E:
 244         case 0x01:
 245             points.append
 246             (
 247                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
 248             );
 249             points.append
 250             (
 251                 generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
 252             );
 253             points.append
 254             (
 255                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
 256             );
 257         break;
 258
 259         case 0x0D:
 260         case 0x02:
 261             points.append
 262             (
 263                 generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
 264             );
 265             points.append
 266             (
 267                 generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
 268             );
 269             points.append
 270             (
 271                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
 272             );
 273         break;
 274
 275         case 0x0C:
 276         case 0x03:
 277         {
 278             Type tp1 =
 279                 generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
 280             Type tp2 =
 281                 generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 282
 283             points.append
 284             (
 285                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
 286             );
 287             points.append(tp1);
 288             points.append(tp2);
 289             points.append(tp2);
 290             points.append
 291             (
 292                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
 293             );
 294             points.append(tp1);
 295         }
 296         break;
 297
 298         case 0x0B:
 299         case 0x04:
 300         {
 301             points.append
 302             (
 303                 generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
 304             );
 305             points.append
 306             (
 307                 generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
 308             );
 309             points.append
 310             (
 311                 generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
 312             );
 313         }
 314         break;
 315
 316         case 0x0A:
 317         case 0x05:
 318         {
 319             Type tp0 =
 320                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
 321             Type tp1 =
 322                 generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 323
 324             points.append(tp0);
 325             points.append(tp1);
 326             points.append
 327             (
 328                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
 329             );
 330             points.append(tp0);
 331             points.append
 332             (
 333                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
 334             );
 335             points.append(tp1);
 336         }
 337         break;
 338
 339         case 0x09:
 340         case 0x06:
 341         {
 342             Type tp0 =
 343                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
 344             Type tp1 =
 345                 generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 346
 347             points.append(tp0);
 348             points.append
 349             (
 350                 generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
 351             );
 352             points.append(tp1);
 353             points.append(tp0);
 354             points.append
 355             (
 356                 generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
 357             );
 358             points.append(tp1);
 359         }
 360         break;
 361
 362         case 0x07:
 363         case 0x08:
 364             points.append
 365             (
 366                 generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
 367             );
 368             points.append
 369             (
 370                 generatePoint(s3,p3,hasSnap3,snapP3,s2,p2,hasSnap2,snapP2)
 371             );
 372             points.append
 373             (
 374                 generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
 375             );
 376         break;
 377     }
 378 }
 379
 380
 381 template<class Type>
 382 Foam::label Foam::isoSurface::generateFaceTriPoints
 383 (
 384     const volScalarField& cVals,
 385     const scalarField& pVals,
 386
 387     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 388     const Field<Type>& pCoords,
 389
 390     const DynamicList<Type>& snappedPoints,
 391     const labelList& snappedCc,
 392     const labelList& snappedPoint,
 393     const label faceI,
 394
 395     const scalar neiVal,
 396     const Type& neiPt,
 397     const bool hasNeiSnap,
 398     const Type& neiSnapPt,
 399
 400     DynamicList<Type>& triPoints,
 401     DynamicList<label>& triMeshCells
 402 ) const
 403 {
 404     label own = mesh_.faceOwner()[faceI];
 405
 406     label oldNPoints = triPoints.size();
 407
 408     const face& f = mesh_.faces()[faceI];
 409
 410     forAll(f, fp)
 411     {
 412         label pointI = f[fp];
 413         label nextPointI = f[f.fcIndex(fp)];
 414
 415         generateTriPoints
 416         (
 417             pVals[pointI],
 418             pCoords[pointI],
 419             snappedPoint[pointI] != -1,
 420             (
 421                 snappedPoint[pointI] != -1
 422               ? snappedPoints[snappedPoint[pointI]]
 423               : pTraits<Type>::zero
 424             ),
 425
 426             pVals[nextPointI],
 427             pCoords[nextPointI],
 428             snappedPoint[nextPointI] != -1,
 429             (
 430                 snappedPoint[nextPointI] != -1
 431               ? snappedPoints[snappedPoint[nextPointI]]
 432               : pTraits<Type>::zero
 433             ),
 434
 435             cVals[own],
 436             cCoords[own],
 437             snappedCc[own] != -1,
 438             (
 439                 snappedCc[own] != -1
 440               ? snappedPoints[snappedCc[own]]
 441               : pTraits<Type>::zero
 442             ),
 443
 444             neiVal,
 445             neiPt,
 446             hasNeiSnap,
 447             neiSnapPt,
 448
 449             triPoints
 450         );
 451     }
 452
 453     // Every three triPoints is a triangle
 454     label nTris = (triPoints.size()-oldNPoints)/3;
 455     for (label i = 0; i < nTris; i++)
 456     {
 457         triMeshCells.append(own);
 458     }
 459
 460     return nTris;
 461 }
 462
 463
 464 template<class Type>
 465 void Foam::isoSurface::generateTriPoints
 466 (
 467     const volScalarField& cVals,
 468     const scalarField& pVals,
 469
 470     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 471     const Field<Type>& pCoords,
 472
 473     const DynamicList<Type>& snappedPoints,
 474     const labelList& snappedCc,
 475     const labelList& snappedPoint,
 476
 477     DynamicList<Type>& triPoints,
 478     DynamicList<label>& triMeshCells
 479 ) const
 480 {
 481     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 482     const labelList& own = mesh_.faceOwner();
 483     const labelList& nei = mesh_.faceNeighbour();
 484
 485     if
 486     (
 487         (cVals.size() != mesh_.nCells())
 488      || (pVals.size() != mesh_.nPoints())
 489      || (cCoords.size() != mesh_.nCells())
 490      || (pCoords.size() != mesh_.nPoints())
 491      || (snappedCc.size() != mesh_.nCells())
 492      || (snappedPoint.size() != mesh_.nPoints())
 493     )
 494     {
 495         FatalErrorIn("isoSurface::generateTriPoints(..)")
 496             << "Incorrect size." << endl
 497             << "mesh: nCells:" << mesh_.nCells()
 498             << " points:" << mesh_.nPoints() << endl
 499             << "cVals:" << cVals.size() << endl
 500             << "cCoords:" << cCoords.size() << endl
 501             << "snappedCc:" << snappedCc.size() << endl
 502             << "pVals:" << pVals.size() << endl
 503             << "pCoords:" << pCoords.size() << endl
 504             << "snappedPoint:" << snappedPoint.size() << endl
 505             << abort(FatalError);
 506     }
 507
 508
 509     // Generate triangle points
 510
 511     triPoints.clear();
 512     triMeshCells.clear();
 513
 514     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 515     {
 516         if (faceCutType_[faceI] != NOTCUT)
 517         {
 518             generateFaceTriPoints
 519             (
 520                 cVals,
 521                 pVals,
 522
 523                 cCoords,
 524                 pCoords,
 525
 526                 snappedPoints,
 527                 snappedCc,
 528                 snappedPoint,
 529                 faceI,
 530
 531                 cVals[nei[faceI]],
 532                 cCoords[nei[faceI]],
 533                 snappedCc[nei[faceI]] != -1,
 534                 (
 535                     snappedCc[nei[faceI]] != -1
 536                   ? snappedPoints[snappedCc[nei[faceI]]]
 537                   : pTraits<Type>::zero
 538                 ),
 539
 540                 triPoints,
 541                 triMeshCells
 542             );
 543         }
 544     }
 545
 546
 547     // Determine neighbouring snap status
 548     boolList neiSnapped(mesh_.nFaces()-mesh_.nInternalFaces(), false);
 549     List<Type> neiSnappedPoint(neiSnapped.size(), pTraits<Type>::zero);
 550     forAll(patches, patchI)
 551     {
 552         const polyPatch& pp = patches[patchI];
 553
 554         if (pp.coupled())
 555         {
 556             label faceI = pp.start();
 557             forAll(pp, i)
 558             {
 559                 label bFaceI = faceI-mesh_.nInternalFaces();
 560                 label snappedIndex = snappedCc[own[faceI]];
 561
 562                 if (snappedIndex != -1)
 563                 {
 564                     neiSnapped[bFaceI] = true;
 565                     neiSnappedPoint[bFaceI] = snappedPoints[snappedIndex];
 566                 }
 567                 faceI++;
 568             }
 569         }
 570     }
 571     syncTools::swapBoundaryFaceList(mesh_, neiSnapped);
 572     syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint);
 573
 574
 575
 576     forAll(patches, patchI)
 577     {
 578         const polyPatch& pp = patches[patchI];
 579
 580         if (isA<processorPolyPatch>(pp))
 581         {
 582             const processorPolyPatch& cpp =
 583                 refCast<const processorPolyPatch>(pp);
 584
 585             PackedBoolList isCollocated(collocatedFaces(cpp));
 586
 587             forAll(isCollocated, i)
 588             {
 589                 label faceI = pp.start()+i;
 590
 591                 if (faceCutType_[faceI] != NOTCUT)
 592                 {
 593                     if (isCollocated[i])
 594                     {
 595                         generateFaceTriPoints
 596                         (
 597                             cVals,
 598                             pVals,
 599
 600                             cCoords,
 601                             pCoords,
 602
 603                             snappedPoints,
 604                             snappedCc,
 605                             snappedPoint,
 606                             faceI,
 607
 608                             cVals.boundaryField()[patchI][i],
 609                             cCoords.boundaryField()[patchI][i],
 610                             neiSnapped[faceI-mesh_.nInternalFaces()],
 611                             neiSnappedPoint[faceI-mesh_.nInternalFaces()],
 612
 613                             triPoints,
 614                             triMeshCells
 615                         );
 616                     }
 617                     else
 618                     {
 619                         generateFaceTriPoints
 620                         (
 621                             cVals,
 622                             pVals,
 623
 624                             cCoords,
 625                             pCoords,
 626
 627                             snappedPoints,
 628                             snappedCc,
 629                             snappedPoint,
 630                             faceI,
 631
 632                             cVals.boundaryField()[patchI][i],
 633                             cCoords.boundaryField()[patchI][i],
 634                             false,
 635                             pTraits<Type>::zero,
 636
 637                             triPoints,
 638                             triMeshCells
 639                         );
 640                     }
 641                 }
 642             }
 643         }
 644         else
 645         {
 646             label faceI = pp.start();
 647
 648             forAll(pp, i)
 649             {
 650                 if (faceCutType_[faceI] != NOTCUT)
 651                 {
 652                     generateFaceTriPoints
 653                     (
 654                         cVals,
 655                         pVals,
 656
 657                         cCoords,
 658                         pCoords,
 659
 660                         snappedPoints,
 661                         snappedCc,
 662                         snappedPoint,
 663                         faceI,
 664
 665                         cVals.boundaryField()[patchI][i],
 666                         cCoords.boundaryField()[patchI][i],
 667                         false,  // fc not snapped
 668                         pTraits<Type>::zero,
 669
 670                         triPoints,
 671                         triMeshCells
 672                     );
 673                 }
 674                 faceI++;
 675             }
 676         }
 677     }
 678
 679     triPoints.shrink();
 680     triMeshCells.shrink();
 681 }
 682
 683
 684 //template <class Type>
 685 //Foam::tmp<Foam::Field<Type> >
 686 //Foam::isoSurface::sample(const Field<Type>& vField) const
 687 //{
 688 //    return tmp<Field<Type> >(new Field<Type>(vField, meshCells()));
 689 //}
 690
 691
 692 template <class Type>
 693 Foam::tmp<Foam::Field<Type> >
 694 Foam::isoSurface::interpolate
 695 (
 696     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 697     const Field<Type>& pCoords
 698 ) const
 699 {
 700     // Recalculate boundary values
 701     tmp<SlicedGeometricField
 702     <
 703         Type,
 704         fvPatchField,
 705         slicedFvPatchField,
 706         volMesh
 707     > > c2(adaptPatchFields(cCoords));
 708
 709
 710     DynamicList<Type> triPoints(nCutCells_);
 711     DynamicList<label> triMeshCells(nCutCells_);
 712
 713     // Dummy snap data
 714     DynamicList<Type> snappedPoints;
 715     labelList snappedCc(mesh_.nCells(), -1);
 716     labelList snappedPoint(mesh_.nPoints(), -1);
 717
 718     generateTriPoints
 719     (
 720         cValsPtr_(),
 721         pVals_,
 722
 723         c2(),
 724         pCoords,
 725
 726         snappedPoints,
 727         snappedCc,
 728         snappedPoint,
 729
 730         triPoints,
 731         triMeshCells
 732     );
 733
 734
 735     // One value per point
 736     tmp<Field<Type> > tvalues
 737     (
 738         new Field<Type>(points().size(), pTraits<Type>::zero)
 739     );
 740     Field<Type>& values = tvalues();
 741     labelList nValues(values.size(), 0);
 742
 743     forAll(triPoints, i)
 744     {
 745         label mergedPointI = triPointMergeMap_[i];
 746
 747         if (mergedPointI >= 0)
 748         {
 749             values[mergedPointI] += triPoints[i];
 750             nValues[mergedPointI]++;
 751         }
 752     }
 753
 754     if (debug)
 755     {
 756         Pout<< "nValues:" << values.size() << endl;
 757         label nMult = 0;
 758         forAll(nValues, i)
 759         {
 760             if (nValues[i] == 0)
 761             {
 762                 FatalErrorIn("isoSurface::interpolate(..)")
 763                     << "point:" << i << " nValues:" << nValues[i]
 764                     << abort(FatalError);
 765             }
 766             else if (nValues[i] > 1)
 767             {
 768                 nMult++;
 769             }
 770         }
 771         Pout<< "Of which mult:" << nMult << endl;
 772     }
 773
 774     forAll(values, i)
 775     {
 776         values[i] /= scalar(nValues[i]);
 777     }
 778
 779     return tvalues;
 780 }
 781
 782
 783 // ************************************************************************* //