src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.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 "sampledTriSurfaceMesh.H"
  27 #include "meshSearch.H"
  28 #include "Tuple2.H"
  29 #include "globalIndex.H"
  30 #include "treeDataPolyMeshCell.H"
  31 #include "treeDataFace.H"
  32 #include "meshTools.H"
  33
  34 #include "addToRunTimeSelectionTable.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(sampledTriSurfaceMesh, 0);
  41     addToRunTimeSelectionTable
  42     (
  43         sampledSurface,
  44         sampledTriSurfaceMesh,
  45         word
  46     );
  47
  48     template<>
  49     const char* NamedEnum<sampledTriSurfaceMesh::samplingSource, 2>::names[] =
  50     {
  51         "cells",
  52         "boundaryFaces"
  53     };
  54
  55     const NamedEnum<sampledTriSurfaceMesh::samplingSource, 2>
  56     sampledTriSurfaceMesh::samplingSourceNames_;
  57
  58
  59     //- Private class for finding nearest
  60     //  Comprising:
  61     //  - global index
  62     //  - sqr(distance)
  63     typedef Tuple2<scalar, label> nearInfo;
  64
  65     class nearestEqOp
  66     {
  67
  68     public:
  69
  70         void operator()(nearInfo& x, const nearInfo& y) const
  71         {
  72             if (y.first() < x.first())
  73             {
  74                 x = y;
  75             }
  76         }
  77     };
  78 }
  79
  80
  81 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  82
  83 const Foam::indexedOctree<Foam::treeDataFace>&
  84 Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
  85 {
  86     // Variant of meshSearch::boundaryTree() that only does non-coupled
  87     // boundary faces.
  88
  89     if (!boundaryTreePtr_.valid())
  90     {
  91         // all non-coupled boundary faces (not just walls)
  92         const polyBoundaryMesh& patches = mesh().boundaryMesh();
  93
  94         labelList bndFaces(mesh().nFaces()-mesh().nInternalFaces());
  95         label bndI = 0;
  96         forAll(patches, patchI)
  97         {
  98             const polyPatch& pp = patches[patchI];
  99             if (!pp.coupled())
 100             {
 101                 forAll(pp, i)
 102                 {
 103                     bndFaces[bndI++] = pp.start()+i;
 104                 }
 105             }
 106         }
 107         bndFaces.setSize(bndI);
 108
 109
 110         treeBoundBox overallBb(mesh().points());
 111         Random rndGen(123456);
 112         overallBb = overallBb.extend(rndGen, 1E-4);
 113         overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 114         overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 115
 116         boundaryTreePtr_.reset
 117         (
 118             new indexedOctree<treeDataFace>
 119             (
 120                 treeDataFace    // all information needed to search faces
 121                 (
 122                     false,                      // do not cache bb
 123                     mesh(),
 124                     bndFaces                    // boundary faces only
 125                 ),
 126                 overallBb,                      // overall search domain
 127                 8,                              // maxLevel
 128                 10,                             // leafsize
 129                 3.0                             // duplicity
 130             )
 131         );
 132     }
 133
 134     return boundaryTreePtr_();
 135 }
 136
 137
 138 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 139
 140 Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 141 (
 142     const word& name,
 143     const polyMesh& mesh,
 144     const word& surfaceName,
 145     const samplingSource sampleSource
 146 )
 147 :
 148     sampledSurface(name, mesh),
 149     surface_
 150     (
 151         IOobject
 152         (
 153             surfaceName,
 154             mesh.time().constant(), // instance
 155             "triSurface",           // local
 156             mesh,                   // registry
 157             IOobject::MUST_READ,
 158             IOobject::NO_WRITE,
 159             false
 160         )
 161     ),
 162     sampleSource_(sampleSource),
 163     needsUpdate_(true),
 164     sampleElements_(0),
 165     samplePoints_(0)
 166 {}
 167
 168
 169 Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 170 (
 171     const word& name,
 172     const polyMesh& mesh,
 173     const dictionary& dict
 174 )
 175 :
 176     sampledSurface(name, mesh, dict),
 177     surface_
 178     (
 179         IOobject
 180         (
 181             dict.lookup("surface"),
 182             mesh.time().constant(), // instance
 183             "triSurface",           // local
 184             mesh,                   // registry
 185             IOobject::MUST_READ,
 186             IOobject::NO_WRITE,
 187             false
 188         )
 189     ),
 190     sampleSource_(samplingSourceNames_[dict.lookup("source")]),
 191     needsUpdate_(true),
 192     sampleElements_(0),
 193     samplePoints_(0)
 194 {}
 195
 196
 197 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 198
 199 Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
 200 {}
 201
 202
 203 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 204
 205 bool Foam::sampledTriSurfaceMesh::needsUpdate() const
 206 {
 207     return needsUpdate_;
 208 }
 209
 210
 211 bool Foam::sampledTriSurfaceMesh::expire()
 212 {
 213     // already marked as expired
 214     if (needsUpdate_)
 215     {
 216         return false;
 217     }
 218
 219     sampledSurface::clearGeom();
 220     MeshStorage::clear();
 221
 222     boundaryTreePtr_.clear();
 223     sampleElements_.clear();
 224     samplePoints_.clear();
 225
 226     needsUpdate_ = true;
 227     return true;
 228 }
 229
 230
 231 bool Foam::sampledTriSurfaceMesh::update()
 232 {
 233     if (!needsUpdate_)
 234     {
 235         return false;
 236     }
 237
 238
 239     // Find the cells the triangles of the surface are in.
 240     // Does approximation by looking at the face centres only
 241     const pointField& fc = surface_.faceCentres();
 242
 243     // Mesh search engine, no triangulation of faces.
 244     meshSearch meshSearcher(mesh(), false);
 245
 246
 247     List<nearInfo> nearest(fc.size());
 248
 249     // Global numbering for cells/faces - only used to uniquely identify local
 250     // elements
 251     globalIndex globalCells
 252     (
 253         sampleSource_ == cells
 254       ? mesh().nCells()
 255       : mesh().nFaces()
 256     );
 257
 258     forAll(nearest, i)
 259     {
 260         nearest[i].first() = GREAT;
 261         nearest[i].second() = labelMax;
 262     }
 263
 264     if (sampleSource_ == cells)
 265     {
 266         // Search for nearest cell
 267
 268         const indexedOctree<treeDataPolyMeshCell>& cellTree =
 269             meshSearcher.cellTree();
 270
 271         forAll(fc, triI)
 272         {
 273             pointIndexHit nearInfo = cellTree.findNearest
 274             (
 275                 fc[triI],
 276                 sqr(GREAT)
 277             );
 278             if (nearInfo.hit())
 279             {
 280                 nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
 281                 nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
 282             }
 283         }
 284     }
 285     else
 286     {
 287         // Search for nearest boundaryFace
 288
 289         ////- Search on all (including coupled) boundary faces
 290         //const indexedOctree<treeDataFace>& bTree = meshSearcher.boundaryTree()
 291         //- Search on all non-coupled boundary faces
 292         const indexedOctree<treeDataFace>& bTree = nonCoupledboundaryTree();
 293
 294         forAll(fc, triI)
 295         {
 296             pointIndexHit nearInfo = bTree.findNearest
 297             (
 298                 fc[triI],
 299                 sqr(GREAT)
 300             );
 301             if (nearInfo.hit())
 302             {
 303                 nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
 304                 nearest[triI].second() = globalCells.toGlobal
 305                 (
 306                     bTree.shapes().faceLabels()[nearInfo.index()]
 307                 );
 308             }
 309         }
 310     }
 311
 312
 313     // See which processor has the nearest. Mark and subset
 314     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 315
 316     Pstream::listCombineGather(nearest, nearestEqOp());
 317     Pstream::listCombineScatter(nearest);
 318
 319     labelList cellOrFaceLabels(fc.size(), -1);
 320
 321     label nFound = 0;
 322     forAll(nearest, triI)
 323     {
 324         if (nearest[triI].second() == labelMax)
 325         {
 326             // Not found on any processor. How to map?
 327         }
 328         else if (globalCells.isLocal(nearest[triI].second()))
 329         {
 330             cellOrFaceLabels[triI] = globalCells.toLocal
 331             (
 332                 nearest[triI].second()
 333             );
 334             nFound++;
 335         }
 336     }
 337
 338
 339     if (debug)
 340     {
 341         Pout<< "Local out of faces:" << cellOrFaceLabels.size()
 342             << " keeping:" << nFound << endl;
 343     }
 344
 345     // Now subset the surface. Do not use triSurface::subsetMesh since requires
 346     // original surface to be in compact numbering.
 347
 348     const triSurface& s = surface_;
 349
 350     // Compact to original triangle
 351     labelList faceMap(s.size());
 352     // Compact to original points
 353     labelList pointMap(s.points().size());
 354     // From original point to compact points
 355     labelList reversePointMap(s.points().size(), -1);
 356
 357     {
 358         label newPointI = 0;
 359         label newFaceI = 0;
 360
 361         forAll(s, faceI)
 362         {
 363             if (cellOrFaceLabels[faceI] != -1)
 364             {
 365                 faceMap[newFaceI++] = faceI;
 366
 367                 const triSurface::FaceType& f = s[faceI];
 368                 forAll(f, fp)
 369                 {
 370                     if (reversePointMap[f[fp]] == -1)
 371                     {
 372                         pointMap[newPointI] = f[fp];
 373                         reversePointMap[f[fp]] = newPointI++;
 374                     }
 375                 }
 376             }
 377         }
 378         faceMap.setSize(newFaceI);
 379         pointMap.setSize(newPointI);
 380     }
 381
 382
 383     // Subset cellOrFaceLabels
 384     cellOrFaceLabels = UIndirectList<label>(cellOrFaceLabels, faceMap)();
 385
 386     // Store any face per point (without using pointFaces())
 387     labelList pointToFace(pointMap.size());
 388
 389     // Create faces and points for subsetted surface
 390     faceList& faces = this->storedFaces();
 391     faces.setSize(faceMap.size());
 392     forAll(faceMap, i)
 393     {
 394         const triFace& f = s[faceMap[i]];
 395         triFace newF
 396         (
 397             reversePointMap[f[0]],
 398             reversePointMap[f[1]],
 399             reversePointMap[f[2]]
 400         );
 401         faces[i] = newF.triFaceFace();
 402
 403         forAll(newF, fp)
 404         {
 405             pointToFace[newF[fp]] = i;
 406         }
 407     }
 408
 409     this->storedPoints() = pointField(s.points(), pointMap);
 410
 411     if (debug)
 412     {
 413         print(Pout);
 414         Pout<< endl;
 415     }
 416
 417
 418
 419     // Collect the samplePoints and sampleElements
 420     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 421
 422     if (sampledSurface::interpolate())
 423     {
 424         samplePoints_.setSize(pointMap.size());
 425         sampleElements_.setSize(pointMap.size(), -1);
 426
 427         if (sampleSource_ == cells)
 428         {
 429             // samplePoints_   : per surface point a location inside the cell
 430             // sampleElements_ : per surface point the cell
 431
 432             forAll(points(), pointI)
 433             {
 434                 const point& pt = points()[pointI];
 435                 label cellI = cellOrFaceLabels[pointToFace[pointI]];
 436                 sampleElements_[pointI] = cellI;
 437
 438                 // Check if point inside cell
 439                 if (mesh().pointInCell(pt, sampleElements_[pointI]))
 440                 {
 441                     samplePoints_[pointI] = pt;
 442                 }
 443                 else
 444                 {
 445                     // Find nearest point on faces of cell
 446                     const cell& cFaces = mesh().cells()[cellI];
 447
 448                     scalar minDistSqr = VGREAT;
 449
 450                     forAll(cFaces, i)
 451                     {
 452                         const face& f = mesh().faces()[cFaces[i]];
 453                         pointHit info = f.nearestPoint(pt, mesh().points());
 454                         if (info.distance() < minDistSqr)
 455                         {
 456                             minDistSqr = info.distance();
 457                             samplePoints_[pointI] = info.rawPoint();
 458                         }
 459                     }
 460                 }
 461             }
 462         }
 463         else
 464         {
 465             // samplePoints_   : per surface point a location on the boundary
 466             // sampleElements_ : per surface point the boundary face containing
 467             //                   the location
 468
 469             forAll(points(), pointI)
 470             {
 471                 const point& pt = points()[pointI];
 472                 label faceI = cellOrFaceLabels[pointToFace[pointI]];
 473                 sampleElements_[pointI] = faceI;
 474                 samplePoints_[pointI] =  mesh().faces()[faceI].nearestPoint
 475                 (
 476                     pt,
 477                     mesh().points()
 478                 ).rawPoint();
 479             }
 480         }
 481     }
 482     else
 483     {
 484         // if sampleSource_ == cells:
 485         //      samplePoints_   : n/a
 486         //      sampleElements_ : per surface triangle the cell
 487         // else:
 488         //      samplePoints_   : n/a
 489         //      sampleElements_ : per surface triangle the boundary face
 490         samplePoints_.clear();
 491         sampleElements_.transfer(cellOrFaceLabels);
 492     }
 493
 494
 495     if (debug)
 496     {
 497         this->clearOut();
 498         OFstream str(mesh().time().path()/"surfaceToMesh.obj");
 499         Info<< "Dumping correspondence from local surface (points or faces)"
 500             << " to mesh (cells or faces) to " << str.name() << endl;
 501         label vertI = 0;
 502
 503         if (sampledSurface::interpolate())
 504         {
 505             if (sampleSource_ == cells)
 506             {
 507                 forAll(samplePoints_, pointI)
 508                 {
 509                     meshTools::writeOBJ(str, points()[pointI]);
 510                     vertI++;
 511
 512                     meshTools::writeOBJ(str, samplePoints_[pointI]);
 513                     vertI++;
 514
 515                     label cellI = sampleElements_[pointI];
 516                     meshTools::writeOBJ(str, mesh().cellCentres()[cellI]);
 517                     vertI++;
 518                     str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
 519                         << nl;
 520                 }
 521             }
 522             else
 523             {
 524                 forAll(samplePoints_, pointI)
 525                 {
 526                     meshTools::writeOBJ(str, points()[pointI]);
 527                     vertI++;
 528
 529                     meshTools::writeOBJ(str, samplePoints_[pointI]);
 530                     vertI++;
 531
 532                     label faceI = sampleElements_[pointI];
 533                     meshTools::writeOBJ(str, mesh().faceCentres()[faceI]);
 534                     vertI++;
 535                     str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
 536                         << nl;
 537                 }
 538             }
 539         }
 540         else
 541         {
 542             if (sampleSource_ == cells)
 543             {
 544                 forAll(sampleElements_, triI)
 545                 {
 546                     meshTools::writeOBJ(str, faceCentres()[triI]);
 547                     vertI++;
 548
 549                     label cellI = sampleElements_[triI];
 550                     meshTools::writeOBJ(str, mesh().cellCentres()[cellI]);
 551                     vertI++;
 552                     str << "l " << vertI-1 << ' ' << vertI << nl;
 553                 }
 554             }
 555             else
 556             {
 557                 forAll(sampleElements_, triI)
 558                 {
 559                     meshTools::writeOBJ(str, faceCentres()[triI]);
 560                     vertI++;
 561
 562                     label faceI = sampleElements_[triI];
 563                     meshTools::writeOBJ(str, mesh().faceCentres()[faceI]);
 564                     vertI++;
 565                     str << "l " << vertI-1 << ' ' << vertI << nl;
 566                 }
 567             }
 568         }
 569     }
 570
 571
 572     needsUpdate_ = false;
 573     return true;
 574 }
 575
 576
 577 Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
 578 (
 579     const volScalarField& vField
 580 ) const
 581 {
 582     return sampleField(vField);
 583 }
 584
 585
 586 Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::sample
 587 (
 588     const volVectorField& vField
 589 ) const
 590 {
 591     return sampleField(vField);
 592 }
 593
 594 Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::sample
 595 (
 596     const volSphericalTensorField& vField
 597 ) const
 598 {
 599     return sampleField(vField);
 600 }
 601
 602
 603 Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::sample
 604 (
 605     const volSymmTensorField& vField
 606 ) const
 607 {
 608     return sampleField(vField);
 609 }
 610
 611
 612 Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::sample
 613 (
 614     const volTensorField& vField
 615 ) const
 616 {
 617     return sampleField(vField);
 618 }
 619
 620
 621 Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::interpolate
 622 (
 623     const interpolation<scalar>& interpolator
 624 ) const
 625 {
 626     return interpolateField(interpolator);
 627 }
 628
 629
 630 Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::interpolate
 631 (
 632     const interpolation<vector>& interpolator
 633 ) const
 634 {
 635     return interpolateField(interpolator);
 636 }
 637
 638 Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::interpolate
 639 (
 640     const interpolation<sphericalTensor>& interpolator
 641 ) const
 642 {
 643     return interpolateField(interpolator);
 644 }
 645
 646
 647 Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::interpolate
 648 (
 649     const interpolation<symmTensor>& interpolator
 650 ) const
 651 {
 652     return interpolateField(interpolator);
 653 }
 654
 655
 656 Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::interpolate
 657 (
 658     const interpolation<tensor>& interpolator
 659 ) const
 660 {
 661     return interpolateField(interpolator);
 662 }
 663
 664
 665 void Foam::sampledTriSurfaceMesh::print(Ostream& os) const
 666 {
 667     os  << "sampledTriSurfaceMesh: " << name() << " :"
 668         << "  surface:" << surface_.objectRegistry::name()
 669         << "  faces:" << faces().size()
 670         << "  points:" << points().size();
 671 }
 672
 673
 674 // ************************************************************************* //