applications/utilities/surface/surfaceToPatch/surfaceToPatch.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     Reads surface and applies surface regioning to a mesh. Uses boundaryMesh
  26     to do the hard work.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "argList.H"
  31 #include "Time.H"
  32 #include "boundaryMesh.H"
  33 #include "polyMesh.H"
  34 #include "faceSet.H"
  35 #include "polyTopoChange.H"
  36 #include "polyModifyFace.H"
  37 #include "globalMeshData.H"
  38
  39 using namespace Foam;
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 // Adds empty patch if not yet there. Returns patchID.
  44 label addPatch(polyMesh& mesh, const word& patchName)
  45 {
  46     label patchI = mesh.boundaryMesh().findPatchID(patchName);
  47
  48     if (patchI == -1)
  49     {
  50         const polyBoundaryMesh& patches = mesh.boundaryMesh();
  51
  52         List<polyPatch*> newPatches(patches.size() + 1);
  53
  54         patchI = 0;
  55
  56         // Copy all old patches
  57         forAll(patches, i)
  58         {
  59             const polyPatch& pp = patches[i];
  60
  61             newPatches[patchI] =
  62                 pp.clone
  63                 (
  64                     patches,
  65                     patchI,
  66                     pp.size(),
  67                     pp.start()
  68                 ).ptr();
  69
  70             patchI++;
  71         }
  72
  73         // Add zero-sized patch
  74         newPatches[patchI] =
  75             new polyPatch
  76             (
  77                 patchName,
  78                 0,
  79                 mesh.nFaces(),
  80                 patchI,
  81                 patches
  82             );
  83
  84         mesh.removeBoundary();
  85         mesh.addPatches(newPatches);
  86
  87         Pout<< "Created patch " << patchName << " at " << patchI << endl;
  88     }
  89     else
  90     {
  91         Pout<< "Reusing patch " << patchName << " at " << patchI << endl;
  92     }
  93
  94     return patchI;
  95 }
  96
  97
  98 // Repatch single face. Return true if patch changed.
  99 bool repatchFace
 100 (
 101     const polyMesh& mesh,
 102     const boundaryMesh& bMesh,
 103     const labelList& nearest,
 104     const labelList& surfToMeshPatch,
 105     const label faceI,
 106     polyTopoChange& meshMod
 107 )
 108 {
 109     bool changed = false;
 110
 111     label bFaceI = faceI - mesh.nInternalFaces();
 112
 113     if (nearest[bFaceI] != -1)
 114     {
 115         // Use boundary mesh one.
 116         label bMeshPatchID = bMesh.whichPatch(nearest[bFaceI]);
 117
 118         label patchID = surfToMeshPatch[bMeshPatchID];
 119
 120         if (patchID != mesh.boundaryMesh().whichPatch(faceI))
 121         {
 122             label own = mesh.faceOwner()[faceI];
 123
 124             label zoneID = mesh.faceZones().whichZone(faceI);
 125
 126             bool zoneFlip = false;
 127
 128             if (zoneID >= 0)
 129             {
 130                 const faceZone& fZone = mesh.faceZones()[zoneID];
 131
 132                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 133             }
 134
 135             meshMod.setAction
 136             (
 137                 polyModifyFace
 138                 (
 139                     mesh.faces()[faceI],// modified face
 140                     faceI,              // label of face being modified
 141                     own,                // owner
 142                     -1,                 // neighbour
 143                     false,              // face flip
 144                     patchID,            // patch for face
 145                     false,              // remove from zone
 146                     zoneID,             // zone for face
 147                     zoneFlip            // face flip in zone
 148                 )
 149             );
 150
 151             changed = true;
 152         }
 153     }
 154     else
 155     {
 156         changed = false;
 157     }
 158     return changed;
 159 }
 160
 161
 162 // Main program:
 163
 164 int main(int argc, char *argv[])
 165 {
 166     argList::addNote
 167     (
 168         "reads surface and applies surface regioning to a mesh"
 169     );
 170
 171     argList::noParallel();
 172     argList::validArgs.append("surfaceFile");
 173     argList::addOption
 174     (
 175         "faceSet",
 176         "name",
 177         "only repatch the faces in specified faceSet"
 178     );
 179     argList::addOption
 180     (
 181         "tol",
 182         "scalar",
 183         "search tolerance as fraction of mesh size (default 1e-3)"
 184     );
 185
 186     #include "setRootCase.H"
 187     #include "createTime.H"
 188     #include "createPolyMesh.H"
 189
 190     const fileName surfName = args[1];
 191
 192     Info<< "Reading surface from " << surfName << " ..." << endl;
 193
 194     word setName;
 195     const bool readSet = args.optionReadIfPresent("faceSet", setName);
 196
 197     if (readSet)
 198     {
 199         Info<< "Repatching only the faces in faceSet " << setName
 200             << " according to nearest surface triangle ..." << endl;
 201     }
 202     else
 203     {
 204         Info<< "Patching all boundary faces according to nearest surface"
 205             << " triangle ..." << endl;
 206     }
 207
 208     const scalar searchTol = args.optionLookupOrDefault("tol", 1e-3);
 209
 210     // Get search box. Anything not within this box will not be considered.
 211     const boundBox& meshBb = mesh.bounds();
 212     const vector searchSpan = searchTol * meshBb.span();
 213
 214     Info<< "All boundary faces further away than " << searchTol
 215         << " of mesh bounding box " << meshBb
 216         << " will keep their patch label ..." << endl;
 217
 218
 219     Info<< "Before patching:" << nl
 220         << "    patch\tsize" << endl;
 221
 222     forAll(mesh.boundaryMesh(), patchI)
 223     {
 224         Info<< "    " << mesh.boundaryMesh()[patchI].name() << '\t'
 225             << mesh.boundaryMesh()[patchI].size() << nl;
 226     }
 227     Info<< endl;
 228
 229
 230     boundaryMesh bMesh;
 231
 232     // Load in the surface.
 233     bMesh.readTriSurface(surfName);
 234
 235     // Add all the boundaryMesh patches to the mesh.
 236     const PtrList<boundaryPatch>& bPatches = bMesh.patches();
 237
 238     // Map from surface patch ( = boundaryMesh patch) to polyMesh patch
 239     labelList patchMap(bPatches.size());
 240
 241     forAll(bPatches, i)
 242     {
 243         patchMap[i] = addPatch(mesh, bPatches[i].name());
 244     }
 245
 246     // Obtain nearest face in bMesh for each boundary face in mesh that
 247     // is within search span.
 248     // Note: should only determine for faceSet if working with that.
 249     labelList nearest(bMesh.getNearest(mesh, searchSpan));
 250
 251     {
 252         // Dump unmatched faces to faceSet for debugging.
 253         faceSet unmatchedFaces(mesh, "unmatchedFaces", nearest.size()/100);
 254
 255         forAll(nearest, bFaceI)
 256         {
 257             if (nearest[bFaceI] == -1)
 258             {
 259                 unmatchedFaces.insert(mesh.nInternalFaces() + bFaceI);
 260             }
 261         }
 262
 263         Pout<< "Writing all " << unmatchedFaces.size()
 264             << " unmatched faces to faceSet "
 265             << unmatchedFaces.name()
 266             << endl;
 267
 268         unmatchedFaces.write();
 269     }
 270
 271
 272     polyTopoChange meshMod(mesh);
 273
 274     label nChanged = 0;
 275
 276     if (readSet)
 277     {
 278         faceSet faceLabels(mesh, setName);
 279         Info<< "Read " << faceLabels.size() << " faces to repatch ..." << endl;
 280
 281         forAllConstIter(faceSet, faceLabels, iter)
 282         {
 283             label faceI = iter.key();
 284
 285             if (repatchFace(mesh, bMesh, nearest, patchMap, faceI, meshMod))
 286             {
 287                 nChanged++;
 288             }
 289         }
 290     }
 291     else
 292     {
 293         forAll(nearest, bFaceI)
 294         {
 295             label faceI = mesh.nInternalFaces() + bFaceI;
 296
 297             if (repatchFace(mesh, bMesh, nearest, patchMap, faceI, meshMod))
 298             {
 299                 nChanged++;
 300             }
 301         }
 302     }
 303
 304     Pout<< "Changed " << nChanged << " boundary faces." << nl << endl;
 305
 306     if (nChanged > 0)
 307     {
 308         meshMod.changeMesh(mesh, false);
 309
 310         Info<< "After patching:" << nl
 311             << "    patch\tsize" << endl;
 312
 313         forAll(mesh.boundaryMesh(), patchI)
 314         {
 315             Info<< "    " << mesh.boundaryMesh()[patchI].name() << '\t'
 316                 << mesh.boundaryMesh()[patchI].size() << endl;
 317         }
 318         Info<< endl;
 319
 320
 321         runTime++;
 322
 323         // Write resulting mesh
 324         Info<< "Writing modified mesh to time " << runTime.value() << endl;
 325         mesh.write();
 326     }
 327
 328
 329     Info<< "End\n" << endl;
 330
 331     return 0;
 332 }
 333
 334
 335 // ************************************************************************* //