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