applications/utilities/mesh/manipulation/autoPatch/autoPatch.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     Divides external faces into patches based on (user supplied) feature
  27     angle.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "argList.H"
  32 #include "polyMesh.H"
  33 #include "Time.H"
  34 #include "boundaryMesh.H"
  35 #include "repatchPolyTopoChanger.H"
  36 #include "mathematicalConstants.H"
  37 #include "OFstream.H"
  38 #include "ListOps.H"
  39
  40 using namespace Foam;
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 // Get all feature edges.
  45 void collectFeatureEdges(const boundaryMesh& bMesh, labelList& markedEdges)
  46 {
  47     markedEdges.setSize(bMesh.mesh().nEdges());
  48
  49     label markedI = 0;
  50
  51     forAll(bMesh.featureSegments(), i)
  52     {
  53         const labelList& segment = bMesh.featureSegments()[i];
  54
  55         forAll(segment, j)
  56         {
  57             label featEdgeI = segment[j];
  58
  59             label meshEdgeI = bMesh.featureToEdge()[featEdgeI];
  60
  61             markedEdges[markedI++] = meshEdgeI;
  62         }
  63     }
  64     markedEdges.setSize(markedI);
  65 }
  66
  67
  68 // Main program:
  69
  70 int main(int argc, char *argv[])
  71 {
  72     argList::noParallel();
  73     argList::validArgs.append("feature angle[0-180]");
  74     argList::validOptions.insert("overwrite", "");
  75
  76 #   include "setRootCase.H"
  77 #   include "createTime.H"
  78     runTime.functionObjects().off();
  79 #   include "createPolyMesh.H"
  80     const word oldInstance = mesh.pointsInstance();
  81
  82     Info<< "Mesh read in = "
  83         << runTime.cpuTimeIncrement()
  84         << " s\n" << endl << endl;
  85
  86
  87     //
  88     // Use boundaryMesh to reuse all the featureEdge stuff in there.
  89     //
  90
  91     boundaryMesh bMesh;
  92
  93     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
  94     bool overwrite = args.optionFound("overwrite");
  95
  96     scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
  97
  98     Info<< "Feature:" << featureAngle << endl
  99         << "minCos :" << minCos << endl
 100         << endl;
 101
 102     bMesh.read(mesh);
 103
 104     // Set feature angle (calculate feature edges)
 105     bMesh.setFeatureEdges(minCos);
 106
 107     // Collect all feature edges as edge labels
 108     labelList markedEdges;
 109
 110     collectFeatureEdges(bMesh, markedEdges);
 111
 112
 113
 114     // (new) patch ID for every face in mesh.
 115     labelList patchIDs(bMesh.mesh().size(), -1);
 116
 117     //
 118     // Fill patchIDs with values for every face by floodfilling without
 119     // crossing feature edge.
 120     //
 121
 122     // Current patch number.
 123     label newPatchI = bMesh.patches().size();
 124
 125     label suffix = 0;
 126
 127     while (true)
 128     {
 129         // Find first unset face.
 130         label unsetFaceI = findIndex(patchIDs, -1);
 131
 132         if (unsetFaceI == -1)
 133         {
 134             // All faces have patchID set. Exit.
 135             break;
 136         }
 137
 138         // Found unset face. Create patch for it.
 139         word patchName;
 140         do
 141         {
 142             patchName = "auto" + name(suffix++);
 143         }
 144         while (bMesh.findPatchID(patchName) != -1);
 145
 146         bMesh.addPatch(patchName);
 147
 148         bMesh.changePatchType(patchName, "patch");
 149
 150
 151         // Fill visited with all faces reachable from unsetFaceI.
 152         boolList visited(bMesh.mesh().size());
 153
 154         bMesh.markFaces(markedEdges, unsetFaceI, visited);
 155
 156
 157         // Assign all visited faces to current patch
 158         label nVisited = 0;
 159
 160         forAll(visited, faceI)
 161         {
 162             if (visited[faceI])
 163             {
 164                 nVisited++;
 165
 166                 patchIDs[faceI] = newPatchI;
 167             }
 168         }
 169
 170         Info<< "Assigned " << nVisited << " faces to patch " << patchName
 171             << endl << endl;
 172
 173         newPatchI++;
 174     }
 175
 176
 177
 178     const PtrList<boundaryPatch>& patches = bMesh.patches();
 179
 180     // Create new list of patches with old ones first
 181     List<polyPatch*> newPatchPtrList(patches.size());
 182
 183     newPatchI = 0;
 184
 185     // Copy old patches
 186     forAll(mesh.boundaryMesh(), patchI)
 187     {
 188         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 189
 190         newPatchPtrList[newPatchI] =
 191             patch.clone
 192             (
 193                 mesh.boundaryMesh(),
 194                 newPatchI,
 195                 patch.size(),
 196                 patch.start()
 197             ).ptr();
 198
 199         newPatchI++;
 200     }
 201
 202     // Add new ones with empty size.
 203     for (label patchI = newPatchI; patchI < patches.size(); patchI++)
 204     {
 205         const boundaryPatch& bp = patches[patchI];
 206
 207         newPatchPtrList[newPatchI] = polyPatch::New
 208         (
 209             polyPatch::typeName,
 210             bp.name(),
 211             0,
 212             mesh.nFaces(),
 213             newPatchI,
 214             mesh.boundaryMesh()
 215         ).ptr();
 216
 217         newPatchI++;
 218     }
 219
 220     if (!overwrite)
 221     {
 222         runTime++;
 223     }
 224
 225
 226     // Change patches
 227     repatchPolyTopoChanger polyMeshRepatcher(mesh);
 228     polyMeshRepatcher.changePatches(newPatchPtrList);
 229
 230
 231     // Change face ordering
 232
 233     // Since bMesh read from mesh there is one to one mapping so we don't
 234     // have to do the geometric stuff.
 235     const labelList& meshFace = bMesh.meshFace();
 236
 237     forAll(patchIDs, faceI)
 238     {
 239         label meshFaceI = meshFace[faceI];
 240
 241         polyMeshRepatcher.changePatchID(meshFaceI, patchIDs[faceI]);
 242     }
 243
 244     polyMeshRepatcher.repatch();
 245
 246     // Write resulting mesh
 247     if (overwrite)
 248     {
 249         mesh.setInstance(oldInstance);
 250     }
 251     mesh.write();
 252
 253
 254     Info<< "End\n" << endl;
 255
 256     return 0;
 257 }
 258
 259
 260 // ************************************************************************* //