applications/utilities/surface/surfaceToPatch/surfaceToPatch.C

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