applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C

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