applications/utilities/surface/surfaceSubset/surfaceSubset.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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     A surface analysis tool which sub-sets the triSurface
  27     to choose only a part of interest. Based on subsetMesh.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "triSurface.H"
  32 #include "argList.H"
  33 #include "OFstream.H"
  34 #include "IFstream.H"
  35 #include "Switch.H"
  36 #include "IOdictionary.H"
  37 #include "boundBox.H"
  38 #include "indexedOctree.H"
  39 #include "octree.H"
  40 #include "treeDataTriSurface.H"
  41 #include "Random.H"
  42
  43 using namespace Foam;
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46 // Main program:
  47
  48 int main(int argc, char *argv[])
  49 {
  50     argList::noParallel();
  51     argList::validArgs.clear();
  52     argList::validArgs.append("surfaceSubsetDict");
  53     argList::validArgs.append("surface file");
  54     argList::validArgs.append("output file");
  55     argList args(argc, argv);
  56
  57     Info<< "Reading dictionary " << args.additionalArgs()[0] << " ..." << endl;
  58     IFstream dictFile(args.additionalArgs()[0]);
  59     dictionary meshSubsetDict(dictFile);
  60
  61     Info<< "Reading surface " << args.additionalArgs()[1] << " ..." << endl;
  62     triSurface surf1(args.additionalArgs()[1]);
  63
  64     Info<< "Original:" << endl;
  65     surf1.writeStats(Info);
  66     Info<< endl;
  67
  68
  69     labelList markedPoints
  70     (
  71         meshSubsetDict.lookup("localPoints")
  72     );
  73
  74     labelList markedEdges
  75     (
  76         meshSubsetDict.lookup("edges")
  77     );
  78
  79     labelList markedFaces
  80     (
  81         meshSubsetDict.lookup("faces")
  82     );
  83
  84     pointField markedZone
  85     (
  86         meshSubsetDict.lookup("zone")
  87     );
  88
  89     if ((markedZone.size() != 0) && (markedZone.size() != 2))
  90     {
  91         FatalErrorIn(args.executable())
  92             << "zone specification should be two points, min and max of "
  93             << "the boundingbox" << endl
  94             << "zone:" << markedZone
  95             << exit(FatalError);
  96     }
  97
  98     Switch addFaceNeighbours
  99     (
 100         meshSubsetDict.lookup("addFaceNeighbours")
 101     );
 102
 103     // Mark the cells for the subset
 104
 105     // Faces to subset
 106     boolList facesToSubset(surf1.size(), false);
 107
 108
 109     //
 110     // pick up faces connected to "localPoints"
 111     //
 112
 113     if (markedPoints.size() > 0)
 114     {
 115         Info << "Found " << markedPoints.size() << " marked point(s)." << endl;
 116
 117         // pick up cells sharing the point
 118
 119         forAll (markedPoints, pointI)
 120         {
 121             if
 122             (
 123                 markedPoints[pointI] < 0
 124              || markedPoints[pointI] >= surf1.nPoints()
 125             )
 126             {
 127                 FatalErrorIn(args.executable())
 128                     << "localPoint label " << markedPoints[pointI]
 129                     << "out of range."
 130                     << " The mesh has got "
 131                     << surf1.nPoints() << " localPoints."
 132                     << exit(FatalError);
 133             }
 134
 135             const labelList& curFaces =
 136                 surf1.pointFaces()[markedPoints[pointI]];
 137
 138             forAll (curFaces, i)
 139             {
 140                 facesToSubset[curFaces[i]] =  true;
 141             }
 142         }
 143     }
 144
 145
 146
 147     //
 148     // pick up faces connected to "edges"
 149     //
 150
 151     if (markedEdges.size() > 0)
 152     {
 153         Info << "Found " << markedEdges.size() << " marked edge(s)." << endl;
 154
 155         // pick up cells sharing the edge
 156
 157         forAll (markedEdges, edgeI)
 158         {
 159             if
 160             (
 161                 markedEdges[edgeI] < 0
 162              || markedEdges[edgeI] >= surf1.nEdges()
 163             )
 164             {
 165                 FatalErrorIn(args.executable())
 166                     << "edge label " << markedEdges[edgeI]
 167                     << "out of range."
 168                     << " The mesh has got "
 169                     << surf1.nEdges() << " edges."
 170                     << exit(FatalError);
 171             }
 172
 173             const labelList& curFaces = surf1.edgeFaces()[markedEdges[edgeI]];
 174
 175             forAll (curFaces, i)
 176             {
 177                 facesToSubset[curFaces[i]] =  true;
 178             }
 179         }
 180     }
 181
 182
 183     //
 184     // pick up faces with centre inside "zone"
 185     //
 186
 187     if (markedZone.size() == 2)
 188     {
 189         const point& min = markedZone[0];
 190         const point& max = markedZone[1];
 191
 192         Info << "Using zone min:" << min << " max:" << max << endl;
 193
 194         forAll(surf1, faceI)
 195         {
 196             const labelledTri& f = surf1[faceI];
 197             const point centre = f.centre(surf1.points());
 198
 199             if
 200             (
 201                 (centre.x() >= min.x())
 202              && (centre.y() >= min.y())
 203              && (centre.z() >= min.z())
 204              && (centre.x() <= max.x())
 205              && (centre.y() <= max.y())
 206              && (centre.z() <= max.z())
 207             )
 208             {
 209                 facesToSubset[faceI] = true;
 210             }
 211         }
 212     }
 213
 214
 215     //
 216     // pick up faces on certain side of surface
 217     //
 218
 219     if (meshSubsetDict.found("surface"))
 220     {
 221         const dictionary& surfDict = meshSubsetDict.subDict("surface");
 222
 223         fileName surfName(surfDict.lookup("name"));
 224
 225         Switch outside(surfDict.lookup("outside"));
 226
 227         if (outside)
 228         {
 229             Info<< "Selecting all triangles with centre outside surface "
 230                 << surfName << endl;
 231         }
 232         else
 233         {
 234             Info<< "Selecting all triangles with centre inside surface "
 235                 << surfName << endl;
 236         }
 237
 238         // Read surface to select on
 239         triSurface selectSurf(surfName);
 240
 241         // bb of surface
 242         treeBoundBox bb(selectSurf.localPoints());
 243
 244         // Radnom number generator
 245         Random rndGen(354543);
 246
 247         // search engine
 248         indexedOctree<treeDataTriSurface> selectTree
 249         (
 250             treeDataTriSurface(selectSurf),
 251             bb.extend(rndGen, 1E-3),    // slightly randomize bb
 252             8,      // maxLevel
 253             10,     // leafsize
 254             3.0     // duplicity
 255         );
 256
 257         // Check if face (centre) is in outside or inside.
 258         forAll(facesToSubset, faceI)
 259         {
 260             if (!facesToSubset[faceI])
 261             {
 262                 const point fc(surf1[faceI].centre(surf1.points()));
 263
 264                 indexedOctree<treeDataTriSurface>::volumeType t =
 265                     selectTree.getVolumeType(fc);
 266
 267                 if (t == indexedOctree<treeDataTriSurface>::INSIDE && !outside)
 268                 {
 269                     facesToSubset[faceI] = true;
 270                 }
 271                 else if
 272                 (
 273                     t == indexedOctree<treeDataTriSurface>::OUTSIDE
 274                  && outside
 275                 )
 276                 {
 277                     facesToSubset[faceI] = true;
 278                 }
 279             }
 280         }
 281     }
 282
 283
 284     //
 285     // pick up specified "faces"
 286     //
 287
 288     // Number of additional faces picked up because of addFaceNeighbours
 289     label nFaceNeighbours = 0;
 290
 291     if (markedFaces.size() > 0)
 292     {
 293         Info << "Found " << markedFaces.size() << " marked face(s)." << endl;
 294
 295         // Check and mark faces to pick up
 296         forAll (markedFaces, faceI)
 297         {
 298             if
 299             (
 300                 markedFaces[faceI] < 0
 301              || markedFaces[faceI] >= surf1.size()
 302             )
 303             {
 304                 FatalErrorIn(args.executable())
 305                     << "Face label " << markedFaces[faceI] << "out of range."
 306                     << " The mesh has got "
 307                     << surf1.size() << " faces."
 308                     << exit(FatalError);
 309             }
 310
 311             // Mark the face
 312             facesToSubset[markedFaces[faceI]] = true;
 313
 314             // mark its neighbours if requested
 315             if (addFaceNeighbours)
 316             {
 317                 const labelList& curFaces =
 318                     surf1.faceFaces()[markedFaces[faceI]];
 319
 320                 forAll (curFaces, i)
 321                 {
 322                     label faceI = curFaces[i];
 323
 324                     if (!facesToSubset[faceI])
 325                     {
 326                         facesToSubset[faceI] =  true;
 327                         nFaceNeighbours++;
 328                     }
 329                 }
 330             }
 331         }
 332     }
 333
 334     if (addFaceNeighbours)
 335     {
 336         Info<< "Added " << nFaceNeighbours
 337             << " faces because of addFaceNeighbours" << endl;
 338     }
 339
 340     // Create subsetted surface
 341     labelList pointMap;
 342     labelList faceMap;
 343     triSurface surf2
 344     (
 345         surf1.subsetMesh(facesToSubset, pointMap, faceMap)
 346     );
 347
 348     Info<< "Subset:" << endl;
 349     surf2.writeStats(Info);
 350     Info << endl;
 351
 352     fileName outFileName(args.additionalArgs()[2]);
 353
 354     Info << "Writing surface to " << outFileName << endl;
 355
 356     surf2.write(outFileName);
 357
 358     return 0;
 359 }
 360
 361
 362 // ************************************************************************* //