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