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) 2004-2010 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
  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     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 "Time.H"
  33 #include "boundaryMesh.H"
  34 #include "repatchPolyTopoChanger.H"
  35 #include "unitConversion.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 #   include "addOverwriteOption.H"
  72     argList::noParallel();
  73     argList::validArgs.append("feature angle[0-180]");
  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     const scalar featureAngle = args.argRead<scalar>(1);
  87     const bool overwrite      = args.optionFound("overwrite");
  88
  89     const scalar minCos = Foam::cos(degToRad(featureAngle));
  90
  91     Info<< "Feature:" << featureAngle << endl
  92         << "minCos :" << minCos << endl
  93         << endl;
  94
  95     //
  96     // Use boundaryMesh to reuse all the featureEdge stuff in there.
  97     //
  98
  99     boundaryMesh bMesh;
 100     bMesh.read(mesh);
 101
 102     // Set feature angle (calculate feature edges)
 103     bMesh.setFeatureEdges(minCos);
 104
 105     // Collect all feature edges as edge labels
 106     labelList markedEdges;
 107
 108     collectFeatureEdges(bMesh, markedEdges);
 109
 110
 111
 112     // (new) patch ID for every face in mesh.
 113     labelList patchIDs(bMesh.mesh().size(), -1);
 114
 115     //
 116     // Fill patchIDs with values for every face by floodfilling without
 117     // crossing feature edge.
 118     //
 119
 120     // Current patch number.
 121     label newPatchI = bMesh.patches().size();
 122
 123     label suffix = 0;
 124
 125     while (true)
 126     {
 127         // Find first unset face.
 128         label unsetFaceI = findIndex(patchIDs, -1);
 129
 130         if (unsetFaceI == -1)
 131         {
 132             // All faces have patchID set. Exit.
 133             break;
 134         }
 135
 136         // Found unset face. Create patch for it.
 137         word patchName;
 138         do
 139         {
 140             patchName = "auto" + name(suffix++);
 141         }
 142         while (bMesh.findPatchID(patchName) != -1);
 143
 144         bMesh.addPatch(patchName);
 145
 146         bMesh.changePatchType(patchName, "patch");
 147
 148
 149         // Fill visited with all faces reachable from unsetFaceI.
 150         boolList visited(bMesh.mesh().size());
 151
 152         bMesh.markFaces(markedEdges, unsetFaceI, visited);
 153
 154
 155         // Assign all visited faces to current patch
 156         label nVisited = 0;
 157
 158         forAll(visited, faceI)
 159         {
 160             if (visited[faceI])
 161             {
 162                 nVisited++;
 163
 164                 patchIDs[faceI] = newPatchI;
 165             }
 166         }
 167
 168         Info<< "Assigned " << nVisited << " faces to patch " << patchName
 169             << endl << endl;
 170
 171         newPatchI++;
 172     }
 173
 174
 175
 176     const PtrList<boundaryPatch>& patches = bMesh.patches();
 177
 178     // Create new list of patches with old ones first
 179     List<polyPatch*> newPatchPtrList(patches.size());
 180
 181     newPatchI = 0;
 182
 183     // Copy old patches
 184     forAll(mesh.boundaryMesh(), patchI)
 185     {
 186         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 187
 188         newPatchPtrList[newPatchI] =
 189             patch.clone
 190             (
 191                 mesh.boundaryMesh(),
 192                 newPatchI,
 193                 patch.size(),
 194                 patch.start()
 195             ).ptr();
 196
 197         newPatchI++;
 198     }
 199
 200     // Add new ones with empty size.
 201     for (label patchI = newPatchI; patchI < patches.size(); patchI++)
 202     {
 203         const boundaryPatch& bp = patches[patchI];
 204
 205         newPatchPtrList[newPatchI] = polyPatch::New
 206         (
 207             polyPatch::typeName,
 208             bp.name(),
 209             0,
 210             mesh.nFaces(),
 211             newPatchI,
 212             mesh.boundaryMesh()
 213         ).ptr();
 214
 215         newPatchI++;
 216     }
 217
 218     if (!overwrite)
 219     {
 220         runTime++;
 221     }
 222
 223
 224     // Change patches
 225     repatchPolyTopoChanger polyMeshRepatcher(mesh);
 226     polyMeshRepatcher.changePatches(newPatchPtrList);
 227
 228
 229     // Change face ordering
 230
 231     // Since bMesh read from mesh there is one to one mapping so we don't
 232     // have to do the geometric stuff.
 233     const labelList& meshFace = bMesh.meshFace();
 234
 235     forAll(patchIDs, faceI)
 236     {
 237         label meshFaceI = meshFace[faceI];
 238
 239         polyMeshRepatcher.changePatchID(meshFaceI, patchIDs[faceI]);
 240     }
 241
 242     polyMeshRepatcher.repatch();
 243
 244     // Write resulting mesh
 245     if (overwrite)
 246     {
 247         mesh.setInstance(oldInstance);
 248     }
 249     mesh.write();
 250
 251
 252     Info<< "End\n" << endl;
 253
 254     return 0;
 255 }
 256
 257
 258 // ************************************************************************* //