applications/utilities/mesh/advanced/selectCells/selectCells.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     Select cells in relation to surface.
  27
  28     Divides cells into three sets:
  29     - cutCells : cells cut by surface or close to surface.
  30     - outside  : cells not in cutCells and reachable from set of
  31       user-defined points (outsidePoints)
  32     - inside   : same but not reachable.
  33
  34     Finally the wanted sets are combined into a cellSet 'selected'. Apart
  35     from straightforward adding the contents there are a few extra rules to
  36     make sure that the surface of the 'outside' of the mesh is singly
  37     connected.
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "argList.H"
  42 #include "Time.H"
  43 #include "polyMesh.H"
  44 #include "IOdictionary.H"
  45 #include "twoDPointCorrector.H"
  46 #include "OFstream.H"
  47 #include "meshTools.H"
  48
  49 #include "triSurface.H"
  50 #include "triSurfaceSearch.H"
  51 #include "meshSearch.H"
  52 #include "cellClassification.H"
  53 #include "cellSet.H"
  54 #include "cellInfo.H"
  55 #include "MeshWave.H"
  56 #include "edgeStats.H"
  57 #include "treeDataTriSurface.H"
  58 #include "indexedOctree.H"
  59
  60 using namespace Foam;
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64 // cellType for cells included/not included in mesh.
  65 static const label MESH = cellClassification::INSIDE;
  66 static const label NONMESH = cellClassification::OUTSIDE;
  67
  68
  69 void writeSet(const cellSet& cells, const string& msg)
  70 {
  71     Info<< "Writing " << msg << " (" << cells.size() << ") to cellSet "
  72         << cells.instance()/cells.local()/cells.name()
  73         << endl << endl;
  74     cells.write();
  75 }
  76
  77
  78 void getType(const labelList& elems, const label type, labelHashSet& set)
  79 {
  80     forAll(elems, i)
  81     {
  82         if (elems[i] == type)
  83         {
  84             set.insert(i);
  85         }
  86     }
  87 }
  88
  89
  90 void cutBySurface
  91 (
  92     const polyMesh& mesh,
  93     const meshSearch& queryMesh,
  94     const triSurfaceSearch& querySurf,
  95
  96     const pointField& outsidePts,
  97     const bool selectCut,
  98     const bool selectInside,
  99     const bool selectOutside,
 100     const scalar nearDist,
 101
 102     cellClassification& cellType
 103 )
 104 {
 105     // Cut with surface and classify as inside/outside/cut
 106     cellType =
 107         cellClassification
 108         (
 109             mesh,
 110             queryMesh,
 111             querySurf,
 112             outsidePts
 113         );
 114
 115     // Get inside/outside/cutCells cellSets.
 116     cellSet inside(mesh, "inside", mesh.nCells()/10);
 117     getType(cellType, cellClassification::INSIDE, inside);
 118     writeSet(inside, "inside cells");
 119
 120     cellSet outside(mesh, "outside", mesh.nCells()/10);
 121     getType(cellType, cellClassification::OUTSIDE, outside);
 122     writeSet(outside, "outside cells");
 123
 124     cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 125     getType(cellType, cellClassification::CUT, cutCells);
 126     writeSet(cutCells, "cells cut by surface");
 127
 128
 129     // Change cellType to reflect selected part of mesh. Use
 130     // MESH to denote selected part, NONMESH for all
 131     // other cells.
 132     // Is a bit of a hack but allows us to reuse all the functionality
 133     // in cellClassification.
 134
 135     forAll(cellType, cellI)
 136     {
 137         label cType = cellType[cellI];
 138
 139         if (cType == cellClassification::CUT)
 140         {
 141             if (selectCut)
 142             {
 143                 cellType[cellI] = MESH;
 144             }
 145             else
 146             {
 147                 cellType[cellI] = NONMESH;
 148             }
 149         }
 150         else if (cType == cellClassification::INSIDE)
 151         {
 152             if (selectInside)
 153             {
 154                 cellType[cellI] = MESH;
 155             }
 156             else
 157             {
 158                 cellType[cellI] = NONMESH;
 159             }
 160         }
 161         else if (cType == cellClassification::OUTSIDE)
 162         {
 163             if (selectOutside)
 164             {
 165                 cellType[cellI] = MESH;
 166             }
 167             else
 168             {
 169                 cellType[cellI] = NONMESH;
 170             }
 171         }
 172         else
 173         {
 174             FatalErrorIn("cutBySurface")
 175                 << "Multiple mesh regions in original mesh" << endl
 176                 << "Please use splitMeshRegions to separate these"
 177                 << exit(FatalError);
 178         }
 179     }
 180
 181
 182     if (nearDist > 0)
 183     {
 184         Info<< "Removing cells with points closer than " << nearDist
 185             << " to the surface ..." << nl << endl;
 186
 187         const pointField& pts = mesh.points();
 188         const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
 189
 190         label nRemoved = 0;
 191
 192         forAll(pts, pointI)
 193         {
 194             const point& pt = pts[pointI];
 195
 196             pointIndexHit hitInfo = tree.findNearest(pt, sqr(nearDist));
 197
 198             if (hitInfo.hit())
 199             {
 200                 const labelList& pCells = mesh.pointCells()[pointI];
 201
 202                 forAll(pCells, i)
 203                 {
 204                     if (cellType[pCells[i]] != NONMESH)
 205                     {
 206                         cellType[pCells[i]] = NONMESH;
 207                         nRemoved++;
 208                     }
 209                 }
 210             }
 211         }
 212
 213 //        tmp<pointField> tnearest = querySurf.calcNearest(pts);
 214 //        const pointField& nearest = tnearest();
 215 //
 216 //        label nRemoved = 0;
 217 //
 218 //        forAll(nearest, pointI)
 219 //        {
 220 //            if (mag(nearest[pointI] - pts[pointI]) < nearDist)
 221 //            {
 222 //                const labelList& pCells = mesh.pointCells()[pointI];
 223 //
 224 //                forAll(pCells, i)
 225 //                {
 226 //                    if (cellType[pCells[i]] != NONMESH)
 227 //                    {
 228 //                        cellType[pCells[i]] = NONMESH;
 229 //                        nRemoved++;
 230 //                    }
 231 //                }
 232 //            }
 233 //        }
 234
 235         Info<< "Removed " << nRemoved << " cells since too close to surface"
 236             << nl << endl;
 237     }
 238 }
 239
 240
 241
 242 // We're meshing the outside. Subset the currently selected mesh cells with the
 243 // ones reachable from the outsidepoints.
 244 label selectOutsideCells
 245 (
 246     const polyMesh& mesh,
 247     const meshSearch& queryMesh,
 248     const pointField& outsidePts,
 249     cellClassification& cellType
 250 )
 251 {
 252     //
 253     // Check all outsidePts and for all of them inside a mesh cell
 254     // collect the faces to start walking from
 255     //
 256
 257     // Outside faces
 258     labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
 259     DynamicList<label> outsideFaces(outsideFacesMap.size());
 260     // CellInfo on outside faces
 261     DynamicList<cellInfo> outsideFacesInfo(outsideFacesMap.size());
 262
 263     // cellInfo for mesh cell
 264     const cellInfo meshInfo(MESH);
 265
 266     forAll(outsidePts, outsidePtI)
 267     {
 268         // Find cell containing point. Linear search.
 269         label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
 270
 271         if (cellType[cellI] == MESH)
 272         {
 273             Info<< "Marking cell " << cellI << " containing outside point "
 274                 << outsidePts[outsidePtI] << " with type " << cellType[cellI]
 275                 << " ..." << endl;
 276
 277             //
 278             // Mark this cell and its faces to start walking from
 279             //
 280
 281             // Mark faces of cellI
 282             const labelList& cFaces = mesh.cells()[cellI];
 283             forAll(cFaces, i)
 284             {
 285                 label faceI = cFaces[i];
 286
 287                 if (outsideFacesMap.insert(faceI))
 288                 {
 289                     outsideFaces.append(faceI);
 290                     outsideFacesInfo.append(meshInfo);
 291                 }
 292             }
 293         }
 294     }
 295
 296     // Floodfill starting from outsideFaces (of type meshInfo)
 297     MeshWave<cellInfo> regionCalc
 298     (
 299         mesh,
 300         outsideFaces.shrink(),
 301         outsideFacesInfo.shrink(),
 302         mesh.nCells()  // max iterations
 303     );
 304
 305     // Now regionCalc should hold info on cells that are reachable from
 306     // changedFaces. Use these to subset cellType
 307     const List<cellInfo>& allCellInfo = regionCalc.allCellInfo();
 308
 309     label nChanged = 0;
 310
 311     forAll(allCellInfo, cellI)
 312     {
 313         if (cellType[cellI] == MESH)
 314         {
 315             // Original cell was selected for meshing. Check if cell was
 316             // reached from outsidePoints
 317             if (allCellInfo[cellI].type() != MESH)
 318             {
 319                 cellType[cellI] = NONMESH;
 320                 nChanged++;
 321             }
 322         }
 323     }
 324
 325     return nChanged;
 326 }
 327
 328
 329 // Main program:
 330
 331 int main(int argc, char *argv[])
 332 {
 333     argList::noParallel();
 334
 335 #   include "setRootCase.H"
 336 #   include "createTime.H"
 337 #   include "createPolyMesh.H"
 338
 339     // Mesh edge statistics calculator
 340     edgeStats edgeCalc(mesh);
 341
 342
 343     IOdictionary refineDict
 344     (
 345         IOobject
 346         (
 347             "selectCellsDict",
 348             runTime.system(),
 349             mesh,
 350             IOobject::MUST_READ,
 351             IOobject::NO_WRITE
 352         )
 353     );
 354
 355     fileName surfName(refineDict.lookup("surface"));
 356     pointField outsidePts(refineDict.lookup("outsidePoints"));
 357     bool useSurface(readBool(refineDict.lookup("useSurface")));
 358     bool selectCut(readBool(refineDict.lookup("selectCut")));
 359     bool selectInside(readBool(refineDict.lookup("selectInside")));
 360     bool selectOutside(readBool(refineDict.lookup("selectOutside")));
 361     scalar nearDist(readScalar(refineDict.lookup("nearDistance")));
 362
 363
 364     if (useSurface)
 365     {
 366         Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
 367             << "    cells cut by surface            : " << selectCut << nl
 368             << "    cells inside of surface         : " << selectInside << nl
 369             << "    cells outside of surface        : " << selectOutside << nl
 370             << "    cells with points further than  : " << nearDist << nl
 371             << endl;
 372     }
 373     else
 374     {
 375         Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
 376             << "    cells reachable from outsidePoints:" << selectOutside << nl
 377             << endl;
 378     }
 379
 380     // Print edge stats on original mesh.
 381     (void)edgeCalc.minLen(Info);
 382
 383     // Search engine on mesh. Face decomposition since faces might be warped.
 384     meshSearch queryMesh(mesh, true);
 385
 386     // Check all 'outside' points
 387     forAll(outsidePts, outsideI)
 388     {
 389         const point& outsidePoint = outsidePts[outsideI];
 390
 391         if (queryMesh.findCell(outsidePoint, -1, false) == -1)
 392         {
 393             FatalErrorIn(args.executable())
 394                 << "outsidePoint " << outsidePoint
 395                 << " is not inside any cell"
 396                 << exit(FatalError);
 397         }
 398     }
 399
 400     // Cell status (compared to surface if provided): inside/outside/cut.
 401     // Start off from everything selected and cut later.
 402     cellClassification cellType
 403     (
 404         mesh,
 405         labelList
 406         (
 407             mesh.nCells(),
 408             cellClassification::MESH
 409         )
 410     );
 411
 412
 413     // Surface
 414     autoPtr<triSurface> surf(NULL);
 415     // Search engine on surface.
 416     autoPtr<triSurfaceSearch> querySurf(NULL);
 417
 418     if (useSurface)
 419     {
 420         surf.reset(new triSurface(surfName));
 421
 422         // Dump some stats
 423         surf().writeStats(Info);
 424
 425         // Search engine on surface.
 426         querySurf.reset(new triSurfaceSearch(surf));
 427
 428         // Set cellType[cellI] according to relation to surface
 429         cutBySurface
 430         (
 431             mesh,
 432             queryMesh,
 433             querySurf,
 434
 435             outsidePts,
 436             selectCut,
 437             selectInside,
 438             selectOutside,
 439             nearDist,
 440
 441             cellType
 442         );
 443     }
 444
 445
 446     // Now 'trim' all the corners from the mesh so meshing/surface extraction
 447     // becomes easier.
 448
 449     label nHanging, nRegionEdges, nRegionPoints, nOutside;
 450
 451     do
 452     {
 453         Info<< "Removing cells which after subsetting would have all points"
 454             << " on outside ..." << nl << endl;
 455
 456         nHanging = cellType.fillHangingCells
 457         (
 458             MESH,       // meshType
 459             NONMESH,    // fill type
 460             mesh.nCells()
 461         );
 462
 463
 464         Info<< "Removing edges connecting cells unconnected by faces ..."
 465             << nl << endl;
 466
 467         nRegionEdges = cellType.fillRegionEdges
 468         (
 469             MESH,       // meshType
 470             NONMESH,    // fill type
 471             mesh.nCells()
 472         );
 473
 474
 475         Info<< "Removing points connecting cells unconnected by faces ..."
 476             << nl << endl;
 477
 478         nRegionPoints = cellType.fillRegionPoints
 479         (
 480             MESH,       // meshType
 481             NONMESH,    // fill type
 482             mesh.nCells()
 483         );
 484
 485         nOutside = 0;
 486         if (selectOutside)
 487         {
 488             // Since we're selecting the cells reachable from outsidePoints
 489             // and the set might have changed, redo the outsideCells
 490             // calculation
 491             nOutside = selectOutsideCells
 492             (
 493                 mesh,
 494                 queryMesh,
 495                 outsidePts,
 496                 cellType
 497             );
 498         }
 499     } while
 500     (
 501         nHanging != 0
 502      || nRegionEdges != 0
 503      || nRegionPoints != 0
 504      || nOutside != 0
 505     );
 506
 507     cellSet selectedCells(mesh, "selected", mesh.nCells()/10);
 508     getType(cellType, MESH, selectedCells);
 509
 510     writeSet(selectedCells, "cells selected for meshing");
 511
 512
 513     Info << "End\n" << endl;
 514
 515     return 0;
 516 }
 517
 518
 519 // ************************************************************************* //