applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.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     Extracts triSurface from a polyMesh. Triangulates all boundary faces.
  26     Region numbers on triangles are the patch numbers of the polyMesh.
  27     Optionally only triangulates named patches.
  28
  29     If run in parallel the processor patches get filtered out by default and
  30     the mesh gets merged. (based on vertex position, not topology, so might go
  31     wrong!).
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "argList.H"
  36 #include "objectRegistry.H"
  37 #include "foamTime.H"
  38 #include "polyMesh.H"
  39 #include "triSurface.H"
  40 #include "triSurfaceTools.H"
  41 #include "processorPolyPatch.H"
  42 #include "ListListOps.H"
  43
  44 using namespace Foam;
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 // Main program:
  49
  50 int main(int argc, char *argv[])
  51 {
  52     argList::validArgs.append("output file");
  53 #   include "addRegionOption.H"
  54     argList::validOptions.insert("excludeProcPatches", "");
  55     argList::validOptions.insert("patches", "(patch0 .. patchN)");
  56
  57 #   include "setRootCase.H"
  58 #   include "createTime.H"
  59
  60     fileName outFileName(runTime.path()/args.additionalArgs()[0]);
  61
  62     Info<< "Extracting triSurface from boundaryMesh ..."
  63         << endl << endl;
  64
  65     Pout<< "Reading mesh from time " << runTime.value() << endl;
  66
  67 #   include "createNamedPolyMesh.H"
  68
  69     bool includeProcPatches =
  70        !(
  71             args.optionFound("excludeProcPatches")
  72          || Pstream::parRun()
  73         );
  74
  75     // Create local surface from:
  76     // - explicitly named patches only (-patches option)
  77     // - all patches (default in sequential mode)
  78     // - all non-processor patches (default in parallel mode)
  79     // - all non-processor patches (sequential mode, -excludeProcPatches option)
  80
  81     // Construct table of patches to include.
  82     const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
  83
  84     labelHashSet includePatches(bMesh.size());
  85
  86     if (args.optionFound("patches"))
  87     {
  88         wordList patchNames(args.optionLookup("patches")());
  89
  90         forAll(patchNames, patchNameI)
  91         {
  92             const word& patchName = patchNames[patchNameI];
  93
  94             label patchI = bMesh.findPatchID(patchName);
  95
  96             if (patchI == -1)
  97             {
  98                 FatalErrorIn(args.executable()) << "No such patch "
  99                     << patchName << endl << "Patches are " << bMesh.names()
 100                     << exit(FatalError);
 101             }
 102             includePatches.insert(patchI);
 103         }
 104     }
 105     else
 106     {
 107         forAll(bMesh, patchI)
 108         {
 109             const polyPatch& patch = bMesh[patchI];
 110
 111             if (includeProcPatches || !isA<processorPolyPatch>(patch))
 112             {
 113                 includePatches.insert(patchI);
 114             }
 115             else
 116             {
 117                 Pout<< patch.name() << " : skipped since processorPatch"
 118                     << endl;
 119             }
 120         }
 121     }
 122
 123     triSurface localSurface
 124     (
 125         triSurfaceTools::triangulate
 126         (
 127             mesh.boundaryMesh(),
 128             includePatches
 129         )
 130     );
 131
 132
 133
 134     if (!Pstream::parRun())
 135     {
 136         Info<< "Writing surface to " << outFileName << endl;
 137
 138         localSurface.write(outFileName);
 139     }
 140     else
 141     {
 142         // Write local surface
 143         fileName localPath = runTime.path()/runTime.caseName() + ".ftr";
 144
 145         Pout<< "Writing local surface to " << localPath << endl;
 146
 147         localSurface.write(localPath);
 148
 149         Info<< endl;
 150
 151
 152         // Gather all points on master
 153         List<pointField> gatheredPoints(Pstream::nProcs());
 154
 155         gatheredPoints[Pstream::myProcNo()] = localSurface.points();
 156
 157         Pstream::gatherList(gatheredPoints);
 158
 159
 160         // Gather all localSurface patches
 161         List<geometricSurfacePatchList> gatheredPatches(Pstream::nProcs());
 162
 163         gatheredPatches[Pstream::myProcNo()] = localSurface.patches();
 164
 165         Pstream::gatherList(gatheredPatches);
 166
 167
 168         // Gather all faces
 169         List<List<labelledTri> > gatheredFaces(Pstream::nProcs());
 170
 171         gatheredFaces[Pstream::myProcNo()] = localSurface;
 172
 173         Pstream::gatherList(gatheredFaces);
 174
 175
 176         if (Pstream::master())
 177         {
 178             // On master combine all points
 179             pointField allPoints =
 180                 ListListOps::combine<pointField>
 181                 (
 182                     gatheredPoints,
 183                     accessOp<pointField>()
 184                 );
 185
 186             // Count number of patches.
 187             label nPatches = 0;
 188
 189             forAll(gatheredPatches, procI)
 190             {
 191                 nPatches += gatheredPatches[procI].size();
 192             }
 193
 194             // Count number of faces.
 195             label nFaces = 0;
 196
 197             forAll(gatheredFaces, procI)
 198             {
 199                 nFaces += gatheredFaces[procI].size();
 200             }
 201
 202
 203
 204             // Loop over all processors and
 205             // - construct mapping from local to global patches
 206             // - relabel faces (both points and regions)
 207
 208             label newPatchI = 0;
 209
 210             // Name to new patchI
 211             HashTable<label> nameToIndex(2*nPatches);
 212
 213             // Storage (oversized) for all patches
 214             geometricSurfacePatchList allPatches(nPatches);
 215
 216             label newFaceI = 0;
 217
 218             // Storage for all faces
 219             List<labelledTri> allFaces(nFaces);
 220
 221             // Offset into allPoints for current processor
 222             label pointOffset = 0;
 223
 224             for (label procI = 0; procI < Pstream::nProcs(); procI++)
 225             {
 226                 Info<< "Processor " << procI << endl
 227                     << "-----------" << endl;
 228
 229                 const geometricSurfacePatchList& patches =
 230                     gatheredPatches[procI];
 231
 232                 // From local patch numbering to global
 233                 labelList localToGlobal(patches.size());
 234
 235                 forAll(patches, patchI)
 236                 {
 237                     const geometricSurfacePatch& sp = patches[patchI];
 238
 239                     if (!nameToIndex.found(sp.name()))
 240                     {
 241                         nameToIndex.insert(sp.name(), newPatchI);
 242
 243                         localToGlobal[patchI] = newPatchI;
 244
 245                         allPatches[newPatchI] = sp;
 246
 247                         newPatchI++;
 248                     }
 249                     else
 250                     {
 251                         localToGlobal[patchI] = nameToIndex[sp.name()];
 252                     }
 253                 }
 254
 255                 Info<< "Local patch to global patch mapping:"
 256                     << endl;
 257
 258                 forAll(patches, patchI)
 259                 {
 260                     Info<< "    name   : " << patches[patchI].name() << endl
 261                         << "    local  : " << patchI << endl
 262                         << "    global : " << localToGlobal[patchI]
 263                         << endl;
 264                 }
 265
 266                 Info<< "Local points added in global points starting at "
 267                     << pointOffset
 268                     << endl;
 269
 270
 271                 // Collect and relabel faces
 272                 const List<labelledTri>& localFaces = gatheredFaces[procI];
 273
 274
 275                 forAll(localFaces, faceI)
 276                 {
 277                     const labelledTri& f = localFaces[faceI];
 278
 279                     allFaces[newFaceI++] =
 280                         labelledTri
 281                         (
 282                             f[0] + pointOffset,
 283                             f[1] + pointOffset,
 284                             f[2] + pointOffset,
 285                             localToGlobal[f.region()]
 286                         );
 287                 }
 288
 289                 pointOffset += gatheredPoints[procI].size();
 290
 291                 Info<< endl;
 292             }
 293             allPatches.setSize(newPatchI);
 294
 295             // We now have allPoints, allFaces and allPatches.
 296             // Construct overall (yet unmerged) surface from these.
 297
 298             triSurface allSurf(allFaces, allPatches, allPoints);
 299
 300             // Cleanup (which does point merge as well
 301             allSurf.cleanup(false);
 302
 303             // Write surface mesh
 304
 305             label slashIndex = runTime.caseName().find_last_of('/');
 306
 307             fileName globalCasePath(runTime.caseName().substr(0, slashIndex));
 308
 309             Info<< "Writing merged surface to " << globalCasePath << endl;
 310
 311             // create database for the sequential run
 312             fileName globalPath
 313             (
 314                 runTime.rootPath()
 315               / globalCasePath
 316               / args.additionalArgs()[0]
 317             );
 318
 319             allSurf.write(globalPath);
 320         }
 321     }
 322
 323     Info << "End\n" << endl;
 324
 325     return 0;
 326 }
 327
 328
 329 // ************************************************************************* //