src/meshTools/surfaceSets/surfaceSets.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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "surfaceSets.H"
  29 #include "polyMesh.H"
  30 #include "triSurface.H"
  31 #include "triSurfaceSearch.H"
  32 #include "pointSet.H"
  33 #include "cellSet.H"
  34 #include "surfaceToCell.H"
  35 #include "cellToPoint.H"
  36 #include "cellToCell.H"
  37 #include "pointToCell.H"
  38 #include "meshSearch.H"
  39 #include "cellClassification.H"
  40
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 //Foam::scalar Foam::surfaceSets::minEdgeLen
  48 //(
  49 //    const primitiveMesh& mesh,
  50 //    const label pointI
  51 //)
  52 //{
  53 //    const edgeList& edges = mesh.edges();
  54 //
  55 //    const pointField& points = mesh.points();
  56 //
  57 //    const labelList& pEdges = mesh.pointEdges()[pointI];
  58 //
  59 //    scalar minLen = GREAT;
  60 //
  61 //    forAll(pEdges, i)
  62 //    {
  63 //        minLen = min(minLen, edges[pEdges[i]].mag(points));
  64 //    }
  65 //    return minLen;
  66 //}
  67 //
  68 //
  69 //// Returns true if cell uses at least one selected point
  70 //bool Foam::surfaceSets::usesPoint
  71 //(
  72 //    const primitiveMesh& mesh,
  73 //    const boolList& selectedPoint,
  74 //    const label cellI
  75 //)
  76 //{
  77 //    const labelList& cFaces = mesh.cells()[cellI];
  78 //
  79 //    forAll(cFaces, cFaceI)
  80 //    {
  81 //        label faceI = cFaces[cFaceI];
  82 //
  83 //        const face& f = mesh.faces()[faceI];
  84 //
  85 //        forAll(f, fp)
  86 //        {
  87 //            if (selectedPoint[f[fp]])
  88 //            {
  89 //                return true;
  90 //            }
  91 //        }
  92 //    }
  93 //    return false;
  94 //}
  95
  96
  97
  98 //// Remove cells in allCells which are connected to other cells in allCells
  99 //// by outside vertices only. Since these outside vertices will be moved onto
 100 //// a surface they might result in flat cells.
 101 //Foam::label Foam::surfaceSets::removeHangingCells
 102 //(
 103 //    const primitiveMesh& mesh,
 104 //    const triSurfaceSearch& querySurf,
 105 //    labelHashSet& internalCells
 106 //)
 107 //{
 108 //    const pointField& points = mesh.points();
 109 //    const cellList& cells = mesh.cells();
 110 //    const faceList& faces = mesh.faces();
 111 //
 112 //    // Determine cells that have all points on the boundary.
 113 //    labelHashSet flatCandidates(getHangingCells(mesh, internalCells));
 114 //
 115 //    // All boundary points will become visible after subsetting and will be
 116 //    // snapped
 117 //    // to surface. Calculate new volume for cells using only these points and
 118 //    // check if it does not become too small.
 119 //
 120 //    // Get points used by flatCandidates.
 121 //    labelHashSet outsidePoints(flatCandidates.size());
 122 //
 123 //    // Snap outside points to surface
 124 //    pointField newPoints(points);
 125 //
 126 //    for
 127 //    (
 128 //        labelHashSet::const_iterator iter = flatCandidates.begin();
 129 //        iter != flatCandidates.end();
 130 //        ++iter
 131 //    )
 132 //    {
 133 //        const cell& cFaces = cells[iter.key()];
 134 //
 135 //        forAll(cFaces, cFaceI)
 136 //        {
 137 //            const face& f = faces[cFaces[cFaceI]];
 138 //
 139 //            forAll(f, fp)
 140 //            {
 141 //                label pointI = f[fp];
 142 //
 143 //                if (outsidePoints.insert(pointI))
 144 //                {
 145 //                    // Calculate new position for this outside point
 146 //                    tmp<pointField> tnearest =
 147 //                        querySurf.calcNearest(pointField(1, points[pointI]));
 148 //                    newPoints[pointI] = tnearest()[0];
 149 //                }
 150 //            }
 151 //        }
 152 //    }
 153 //
 154 //
 155 //    // Calculate new volume for mixed cells
 156 //    label nRemoved = 0;
 157 //    for
 158 //    (
 159 //        labelHashSet::const_iterator iter = flatCandidates.begin();
 160 //        iter != flatCandidates.end();
 161 //        ++iter
 162 //    )
 163 //    {
 164 //        label cellI = iter.key();
 165 //
 166 //        const cell& cll = cells[cellI];
 167 //
 168 //        scalar newVol = cll.mag(newPoints, faces);
 169 //        scalar oldVol = mesh.cellVolumes()[cellI];
 170 //
 171 //        if (newVol < 0.1 * oldVol)
 172 //        {
 173 //            internalCells.erase(cellI);
 174 //            nRemoved++;
 175 //        }
 176 //    }
 177 //
 178 //    return nRemoved;
 179 //}
 180
 181
 182 //// Select all points out of pointSet where the distance to the surface is less
 183 //// than a factor times a local length scale (minimum length of connected edges)
 184 //void Foam::surfaceSets::getNearPoints
 185 //(
 186 //    const primitiveMesh& mesh,
 187 //    const triSurface&,
 188 //    const triSurfaceSearch& querySurf,
 189 //    const scalar edgeFactor,
 190 //    const pointSet& candidateSet,
 191 //    pointSet& nearPointSet
 192 //)
 193 //{
 194 //    if (edgeFactor <= 0)
 195 //    {
 196 //        return;
 197 //    }
 198 //
 199 //    labelList candidates(candidateSet.toc());
 200 //
 201 //    pointField candidatePoints(candidates.size());
 202 //    forAll(candidates, i)
 203 //    {
 204 //        candidatePoints[i] = mesh.points()[candidates[i]];
 205 //    }
 206 //
 207 //    tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
 208 //    const pointField& nearest = tnearest();
 209 //
 210 //    const pointField& points = mesh.points();
 211 //
 212 //    forAll(candidates, i)
 213 //    {
 214 //        label pointI = candidates[i];
 215 //
 216 //        scalar minLen = minEdgeLen(mesh, pointI);
 217 //
 218 //        scalar dist = mag(nearest[i] - points[pointI]);
 219 //
 220 //        if (dist < edgeFactor * minLen)
 221 //        {
 222 //            nearPointSet.insert(pointI);
 223 //        }
 224 //    }
 225 //}
 226
 227
 228 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 229
 230
 231 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 232
 233 void Foam::surfaceSets::getSurfaceSets
 234 (
 235     const polyMesh& mesh,
 236     const fileName&,
 237     const triSurface&,
 238     const triSurfaceSearch& querySurf,
 239     const pointField& outsidePts,
 240
 241     const label nCutLayers,
 242
 243     labelHashSet& inside,
 244     labelHashSet& outside,
 245     labelHashSet& cut
 246 )
 247 {
 248     // Construct search engine on mesh
 249     meshSearch queryMesh(mesh, true);
 250
 251
 252     // Check all 'outside' points
 253     forAll(outsidePts, outsideI)
 254     {
 255         const point& outsidePoint = outsidePts[outsideI];
 256
 257         // Find cell point is in. Linear search.
 258         if (queryMesh.findCell(outsidePoint, -1, false) == -1)
 259         {
 260             FatalErrorIn
 261             (
 262                 "surfaceSets::getSurfaceSets"
 263                 "(const polyMesh&, const fileName&, const triSurface&"
 264                 ", const triSurfaceSearch&, const pointField&"
 265                 ", labelHashSet&, labelHashSet&, labelHashSet&)"
 266             )   << "outsidePoint " << outsidePoint
 267                 << " is not inside any cell"
 268                 << exit(FatalError);
 269         }
 270     }
 271
 272     // Cut faces with surface and classify cells
 273     cellClassification cellType
 274     (
 275         mesh,
 276         queryMesh,
 277         querySurf,
 278         outsidePts
 279     );
 280
 281     if (nCutLayers > 0)
 282     {
 283         // Trim cutCells so they are max nCutLayers away (as seen in point-cell
 284         // walk) from outside cells.
 285         cellType.trimCutCells
 286         (
 287             nCutLayers,
 288             cellClassification::OUTSIDE,
 289             cellClassification::INSIDE
 290         );
 291     }
 292
 293     forAll(cellType, cellI)
 294     {
 295         label cType = cellType[cellI];
 296
 297         if (cType == cellClassification::CUT)
 298         {
 299             cut.insert(cellI);
 300         }
 301         else if (cType == cellClassification::INSIDE)
 302         {
 303             inside.insert(cellI);
 304         }
 305         else if (cType == cellClassification::OUTSIDE)
 306         {
 307             outside.insert(cellI);
 308         }
 309     }
 310 }
 311
 312
 313 Foam::labelHashSet Foam::surfaceSets::getHangingCells
 314 (
 315     const primitiveMesh& mesh,
 316     const labelHashSet& internalCells
 317 )
 318 {
 319     const cellList& cells = mesh.cells();
 320     const faceList& faces = mesh.faces();
 321
 322
 323     // Divide points into
 324     // -referenced by internal only
 325     // -referenced by outside only
 326     // -mixed
 327
 328     List<pointStatus> pointSide(mesh.nPoints(), NOTSET);
 329
 330     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 331     {
 332         if (internalCells.found(cellI))
 333         {
 334             // Inside cell. Mark all vertices seen from this cell.
 335             const labelList& cFaces = cells[cellI];
 336
 337             forAll(cFaces, cFaceI)
 338             {
 339                 const face& f = faces[cFaces[cFaceI]];
 340
 341                 forAll(f, fp)
 342                 {
 343                     label pointI = f[fp];
 344
 345                     if (pointSide[pointI] == NOTSET)
 346                     {
 347                         pointSide[pointI] = INSIDE;
 348                     }
 349                     else if (pointSide[pointI] == OUTSIDE)
 350                     {
 351                         pointSide[pointI] = MIXED;
 352                     }
 353                     else
 354                     {
 355                         // mixed or inside stay same
 356                     }
 357                 }
 358             }
 359         }
 360         else
 361         {
 362             // Outside cell
 363             const labelList& cFaces = cells[cellI];
 364
 365             forAll(cFaces, cFaceI)
 366             {
 367                 const face& f = faces[cFaces[cFaceI]];
 368
 369                 forAll(f, fp)
 370                 {
 371                     label pointI = f[fp];
 372
 373                     if (pointSide[pointI] == NOTSET)
 374                     {
 375                         pointSide[pointI] = OUTSIDE;
 376                     }
 377                     else if (pointSide[pointI] == INSIDE)
 378                     {
 379                         pointSide[pointI] = MIXED;
 380                     }
 381                     else
 382                     {
 383                         // mixed or outside stay same
 384                     }
 385                 }
 386             }
 387         }
 388     }
 389
 390
 391     //OFstream mixedStr("mixed.obj");
 392     //
 393     //forAll(pointSide, pointI)
 394     //{
 395     //    if (pointSide[pointI] == MIXED)
 396     //    {
 397     //        const point& pt = points[pointI];
 398     //
 399     //        mixedStr << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
 400     //            << endl;
 401     //    }
 402     //}
 403
 404
 405     // Determine cells using mixed points only
 406
 407     labelHashSet mixedOnlyCells(internalCells.size());
 408
 409     for
 410     (
 411         labelHashSet::const_iterator iter = internalCells.begin();
 412         iter != internalCells.end();
 413         ++iter
 414     )
 415     {
 416         label cellI = iter.key();
 417
 418         const cell& cFaces = cells[cellI];
 419
 420         label usesMixedOnly = true;
 421
 422         forAll(cFaces, i)
 423         {
 424             const face& f = faces[cFaces[i]];
 425
 426             forAll(f, fp)
 427             {
 428                 if (pointSide[f[fp]] != MIXED)
 429                 {
 430                     usesMixedOnly = false;
 431                     break;
 432                 }
 433             }
 434
 435             if (!usesMixedOnly)
 436             {
 437                 break;
 438             }
 439         }
 440         if (usesMixedOnly)
 441         {
 442             mixedOnlyCells.insert(cellI);
 443         }
 444     }
 445
 446     return mixedOnlyCells;
 447 }
 448
 449
 450 //void Foam::surfaceSets::writeSurfaceSets
 451 //(
 452 //    const polyMesh& mesh,
 453 //    const fileName& surfName,
 454 //    const triSurface& surf,
 455 //    const triSurfaceSearch& querySurf,
 456 //    const pointField& outsidePts,
 457 //    const scalar edgeFactor
 458 //)
 459 //{
 460 //    // Cellsets for inside/outside determination
 461 //    cellSet rawInside(mesh, "rawInside", mesh.nCells()/10);
 462 //    cellSet rawOutside(mesh, "rawOutside", mesh.nCells()/10);
 463 //    cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 464 //
 465 //    // Get inside/outside/cut cells
 466 //    getSurfaceSets
 467 //    (
 468 //        mesh,
 469 //        surfName,
 470 //        surf,
 471 //        querySurf,
 472 //        outsidePts,
 473 //
 474 //        rawInside,
 475 //        rawOutside,
 476 //        cutCells
 477 //    );
 478 //
 479 //
 480 //    Pout<< "rawInside:" << rawInside.size() << endl;
 481 //
 482 //    label nRemoved;
 483 //    do
 484 //    {
 485 //        nRemoved = removeHangingCells(mesh, querySurf, rawInside);
 486 //
 487 //        Pout<< nl
 488 //            << "Removed " << nRemoved
 489 //            << " rawInside cells that have all their points on the outside"
 490 //            << endl;
 491 //    }
 492 //    while (nRemoved != 0);
 493 //
 494 //    Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
 495 //        << rawInside.instance()/rawInside.local()/rawInside.name()
 496 //        << endl << endl;
 497 //    rawInside.write();
 498 //
 499 //
 500 //
 501 //    // Select outside cells
 502 //    surfaceToCell outsideSource
 503 //    (
 504 //        mesh,
 505 //        surfName,
 506 //        surf,
 507 //        querySurf,
 508 //        outsidePts,
 509 //        false,          // includeCut
 510 //        false,          // includeInside
 511 //        true,           // includeOutside
 512 //        -GREAT,         // nearDist
 513 //        -GREAT          // curvature
 514 //    );
 515 //
 516 //    outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
 517 //
 518 //    Pout<< "rawOutside:" << rawOutside.size() << endl;
 519 //
 520 //    do
 521 //    {
 522 //        nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
 523 //
 524 //        Pout<< nl
 525 //            << "Removed " << nRemoved
 526 //            << " rawOutside cells that have all their points on the outside"
 527 //            << endl;
 528 //    }
 529 //    while (nRemoved != 0);
 530 //
 531 //    Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
 532 //        << rawOutside.instance()/rawOutside.local()/rawOutside.name()
 533 //        << endl << endl;
 534 //    rawOutside.write();
 535 //
 536 //
 537 //    // Select cut cells by negating inside and outside set.
 538 //    cutCells.invert(mesh.nCells());
 539 //
 540 //    cellToCell deleteInsideSource(mesh, rawInside.name());
 541 //
 542 //    deleteInsideSource.applyToSet(topoSetSource::DELETE, cutCells);
 543 //    Pout<< "Writing cut cells (" << cutCells.size() << ") to cellSet "
 544 //        << cutCells.instance()/cutCells.local()/cutCells.name()
 545 //        << endl << endl;
 546 //    cutCells.write();
 547 //
 548 //
 549 //    //
 550 //    // Remove cells with points too close to surface.
 551 //    //
 552 //
 553 //
 554 //    // Get all points in cutCells.
 555 //    pointSet cutPoints(mesh, "cutPoints", 4*cutCells.size());
 556 //    cellToPoint cutSource(mesh, "cutCells", cellToPoint::ALL);
 557 //    cutSource.applyToSet(topoSetSource::NEW, cutPoints);
 558 //
 559 //    // Get all points that are too close to surface.
 560 //    pointSet nearPoints(mesh, "nearPoints", cutPoints.size());
 561 //
 562 //    getNearPoints
 563 //    (
 564 //        mesh,
 565 //        surf,
 566 //        querySurf,
 567 //        edgeFactor,
 568 //        cutPoints,
 569 //        nearPoints
 570 //    );
 571 //
 572 //    Pout<< nl
 573 //        << "Selected " << nearPoints.size()
 574 //        << " points that are closer than " << edgeFactor
 575 //        << " times the local minimum lengthscale to the surface"
 576 //        << nl << endl;
 577 //
 578 //
 579 //    // Remove cells that use any of the points in nearPoints
 580 //    // from the inside and outsideCells.
 581 //    nearPoints.write();
 582 //    pointToCell pToCell(mesh, nearPoints.name(), pointToCell::ANY);
 583 //
 584 //
 585 //
 586 //    // Start off from copy of rawInside, rawOutside
 587 //    cellSet inside(mesh, "inside", rawInside);
 588 //    cellSet outside(mesh, "outside", rawOutside);
 589 //
 590 //    pToCell.applyToSet(topoSetSource::DELETE, inside);
 591 //    pToCell.applyToSet(topoSetSource::DELETE, outside);
 592 //
 593 //    Pout<< nl
 594 //        << "Removed " << rawInside.size() - inside.size()
 595 //        << " inside cells that are too close to the surface" << endl;
 596 //
 597 //    Pout<< nl
 598 //        << "Removed " << rawOutside.size() - outside.size()
 599 //        << " inside cells that are too close to the surface" << nl << endl;
 600 //
 601 //
 602 //
 603 //    //
 604 //    // Remove cells with one or no faces on rest of cellSet. Note: Problem is
 605 //    // not these cells an sich but rather that all of their vertices will be
 606 //    // outside vertices and thus projected onto surface. Which might (if they
 607 //    // project onto same surface) result in flattened cells.
 608 //    //
 609 //
 610 //    do
 611 //    {
 612 //        nRemoved = removeHangingCells(mesh, querySurf, inside);
 613 //
 614 //        Pout<< nl
 615 //            << "Removed " << nRemoved
 616 //            << " inside cells that have all their points on the outside"
 617 //            << endl;
 618 //    }
 619 //    while (nRemoved != 0);
 620 //    do
 621 //    {
 622 //        nRemoved = removeHangingCells(mesh, querySurf, outside);
 623 //
 624 //        Pout<< nl
 625 //            << "Removed " << nRemoved
 626 //            << " outside cells that have all their points on the inside"
 627 //            << endl;
 628 //    }
 629 //    while (nRemoved != 0);
 630 //
 631 //
 632 //    //
 633 //    // Write
 634 //    //
 635 //
 636 //
 637 //    Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
 638 //        << inside.instance()/inside.local()/inside.name()
 639 //        << endl << endl;
 640 //
 641 //    inside.write();
 642 //
 643 //    Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
 644 //        << outside.instance()/outside.local()/outside.name()
 645 //        << endl << endl;
 646 //
 647 //    outside.write();
 648 //}
 649
 650
 651 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 652
 653
 654 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 655
 656
 657 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 658
 659
 660 // ************************************************************************* //