applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C

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