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