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