applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.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 Application
  25     mergeOrSplitBaffles
  26
  27 Description
  28     Detects faces that share points (baffles). Either merge them or
  29     duplicate the points.
  30
  31     Notes:
  32     - can only handle pairwise boundary faces. So three faces using
  33       the same points is not handled (is illegal mesh anyway)
  34
  35     - there is no option to only split/merge some baffles.
  36
  37     - surfaces consisting of duplicate faces can be topologically split
  38     if the points on the interior of the surface cannot walk to all the
  39     cells that use them in one go.
  40
  41     - Parallel operation (where duplicate face is perpendicular to a coupled
  42     boundary) is supported but not really tested.
  43     (Note that coupled faces themselves are not seen as duplicate faces)
  44
  45 \*---------------------------------------------------------------------------*/
  46
  47 #include "argList.H"
  48 #include "Time.H"
  49 #include "syncTools.H"
  50 #include "faceSet.H"
  51 #include "pointSet.H"
  52 #include "meshTools.H"
  53 #include "polyTopoChange.H"
  54 #include "polyRemoveFace.H"
  55 #include "polyModifyFace.H"
  56 #include "indirectPrimitivePatch.H"
  57 #include "processorPolyPatch.H"
  58 #include "localPointRegion.H"
  59 #include "duplicatePoints.H"
  60 #include "ReadFields.H"
  61 #include "volFields.H"
  62 #include "surfaceFields.H"
  63
  64 using namespace Foam;
  65
  66 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  67
  68 void insertDuplicateMerge
  69 (
  70     const polyMesh& mesh,
  71     const labelList& duplicates,
  72     polyTopoChange& meshMod
  73 )
  74 {
  75     const faceList& faces = mesh.faces();
  76     const labelList& faceOwner = mesh.faceOwner();
  77     const faceZoneMesh& faceZones = mesh.faceZones();
  78
  79     forAll(duplicates, bFaceI)
  80     {
  81         label otherFaceI = duplicates[bFaceI];
  82
  83         if (otherFaceI != -1 && otherFaceI > bFaceI)
  84         {
  85             // Two duplicate faces. Merge.
  86
  87             label face0 = mesh.nInternalFaces() + bFaceI;
  88             label face1 = mesh.nInternalFaces() + otherFaceI;
  89
  90             label own0 = faceOwner[face0];
  91             label own1 = faceOwner[face1];
  92
  93             if (own0 < own1)
  94             {
  95                 // Use face0 as the new internal face.
  96                 label zoneID = faceZones.whichZone(face0);
  97                 bool zoneFlip = false;
  98
  99                 if (zoneID >= 0)
 100                 {
 101                     const faceZone& fZone = faceZones[zoneID];
 102                     zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
 103                 }
 104
 105                 meshMod.setAction(polyRemoveFace(face1));
 106                 meshMod.setAction
 107                 (
 108                     polyModifyFace
 109                     (
 110                         faces[face0],           // modified face
 111                         face0,                  // label of face being modified
 112                         own0,                   // owner
 113                         own1,                   // neighbour
 114                         false,                  // face flip
 115                         -1,                     // patch for face
 116                         false,                  // remove from zone
 117                         zoneID,                 // zone for face
 118                         zoneFlip                // face flip in zone
 119                     )
 120                 );
 121             }
 122             else
 123             {
 124                 // Use face1 as the new internal face.
 125                 label zoneID = faceZones.whichZone(face1);
 126                 bool zoneFlip = false;
 127
 128                 if (zoneID >= 0)
 129                 {
 130                     const faceZone& fZone = faceZones[zoneID];
 131                     zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
 132                 }
 133
 134                 meshMod.setAction(polyRemoveFace(face0));
 135                 meshMod.setAction
 136                 (
 137                     polyModifyFace
 138                     (
 139                         faces[face1],           // modified face
 140                         face1,                  // label of face being modified
 141                         own1,                   // owner
 142                         own0,                   // neighbour
 143                         false,                  // face flip
 144                         -1,                     // patch for face
 145                         false,                  // remove from zone
 146                         zoneID,                 // zone for face
 147                         zoneFlip                // face flip in zone
 148                     )
 149                 );
 150             }
 151         }
 152     }
 153 }
 154
 155
 156 labelList findBaffles(const polyMesh& mesh, const labelList& boundaryFaces)
 157 {
 158     // Get all duplicate face labels (in boundaryFaces indices!).
 159     labelList duplicates = localPointRegion::findDuplicateFaces
 160     (
 161         mesh,
 162         boundaryFaces
 163     );
 164
 165
 166     // Check that none are on processor patches
 167     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 168
 169     forAll(duplicates, bFaceI)
 170     {
 171         if (duplicates[bFaceI] != -1)
 172         {
 173             label faceI = mesh.nInternalFaces() + bFaceI;
 174             label patchI = patches.whichPatch(faceI);
 175
 176             if (isA<processorPolyPatch>(patches[patchI]))
 177             {
 178                 FatalErrorIn("findBaffles(const polyMesh&, const labelList&)")
 179                     << "Duplicate face " << faceI
 180                     << " is on a processorPolyPatch."
 181                     << "This is not allowed." << nl
 182                     << "Face:" << faceI
 183                     << " is on patch:" << patches[patchI].name()
 184                     << abort(FatalError);
 185             }
 186         }
 187     }
 188
 189
 190     // Write to faceSet for ease of postprocessing.
 191     {
 192         faceSet duplicateSet
 193         (
 194             mesh,
 195             "duplicateFaces",
 196             (mesh.nFaces() - mesh.nInternalFaces())/256
 197         );
 198
 199         forAll(duplicates, bFaceI)
 200         {
 201             label otherFaceI = duplicates[bFaceI];
 202
 203             if (otherFaceI != -1 && otherFaceI > bFaceI)
 204             {
 205                 duplicateSet.insert(mesh.nInternalFaces() + bFaceI);
 206                 duplicateSet.insert(mesh.nInternalFaces() + otherFaceI);
 207             }
 208         }
 209
 210         Pout<< "Writing " << duplicateSet.size()
 211             << " duplicate faces to faceSet " << duplicateSet.objectPath()
 212             << nl << endl;
 213         duplicateSet.write();
 214     }
 215
 216     return duplicates;
 217 }
 218
 219
 220
 221
 222 int main(int argc, char *argv[])
 223 {
 224     argList::addNote
 225     (
 226         "Detect faces that share points (baffles).\n"
 227         "Merge them or duplicate the points."
 228     );
 229
 230     #include "addOverwriteOption.H"
 231     #include "addRegionOption.H"
 232     argList::addBoolOption
 233     (
 234         "detectOnly",
 235         "find baffles only, but do not merge or split them"
 236     );
 237     argList::addBoolOption
 238     (
 239         "split",
 240         "topologically split duplicate surfaces"
 241     );
 242
 243     #include "setRootCase.H"
 244     #include "createTime.H"
 245     runTime.functionObjects().off();
 246     #include "createNamedMesh.H"
 247
 248     const word oldInstance = mesh.pointsInstance();
 249
 250     const bool split      = args.optionFound("split");
 251     const bool overwrite  = args.optionFound("overwrite");
 252     const bool detectOnly = args.optionFound("detectOnly");
 253
 254     // Collect all boundary faces
 255     labelList boundaryFaces(mesh.nFaces() - mesh.nInternalFaces());
 256
 257     forAll(boundaryFaces, i)
 258     {
 259         boundaryFaces[i] = i+mesh.nInternalFaces();
 260     }
 261
 262
 263     if (detectOnly)
 264     {
 265         findBaffles(mesh, boundaryFaces);
 266         return 0;
 267     }
 268
 269
 270     // Read objects in time directory
 271     IOobjectList objects(mesh, runTime.timeName());
 272
 273     // Read vol fields.
 274
 275     PtrList<volScalarField> vsFlds;
 276     ReadFields(mesh, objects, vsFlds);
 277
 278     PtrList<volVectorField> vvFlds;
 279     ReadFields(mesh, objects, vvFlds);
 280
 281     PtrList<volSphericalTensorField> vstFlds;
 282     ReadFields(mesh, objects, vstFlds);
 283
 284     PtrList<volSymmTensorField> vsymtFlds;
 285     ReadFields(mesh, objects, vsymtFlds);
 286
 287     PtrList<volTensorField> vtFlds;
 288     ReadFields(mesh, objects, vtFlds);
 289
 290     // Read surface fields.
 291
 292     PtrList<surfaceScalarField> ssFlds;
 293     ReadFields(mesh, objects, ssFlds);
 294
 295     PtrList<surfaceVectorField> svFlds;
 296     ReadFields(mesh, objects, svFlds);
 297
 298     PtrList<surfaceSphericalTensorField> sstFlds;
 299     ReadFields(mesh, objects, sstFlds);
 300
 301     PtrList<surfaceSymmTensorField> ssymtFlds;
 302     ReadFields(mesh, objects, ssymtFlds);
 303
 304     PtrList<surfaceTensorField> stFlds;
 305     ReadFields(mesh, objects, stFlds);
 306
 307
 308     // Mesh change engine
 309     polyTopoChange meshMod(mesh);
 310
 311
 312     if (split)
 313     {
 314         Pout<< "Topologically splitting duplicate surfaces"
 315             << ", i.e. duplicating points internal to duplicate surfaces."
 316             << nl << endl;
 317
 318         // Analyse which points need to be duplicated
 319         localPointRegion regionSide(mesh);
 320
 321         // Point duplication engine
 322         duplicatePoints pointDuplicator(mesh);
 323
 324         // Insert topo changes
 325         pointDuplicator.setRefinement(regionSide, meshMod);
 326     }
 327     else
 328     {
 329         Pout<< "Merging duplicate faces."
 330             << nl << endl;
 331
 332         // Get all duplicate face labels (in boundaryFaces indices!).
 333         labelList duplicates(findBaffles(mesh, boundaryFaces));
 334
 335         // Merge into internal faces.
 336         insertDuplicateMerge(mesh, duplicates, meshMod);
 337     }
 338
 339     if (!overwrite)
 340     {
 341         runTime++;
 342     }
 343
 344     // Change the mesh. No inflation.
 345     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
 346
 347     // Update fields
 348     mesh.updateMesh(map);
 349
 350     // Move mesh (since morphing does not do this)
 351     if (map().hasMotionPoints())
 352     {
 353         mesh.movePoints(map().preMotionPoints());
 354     }
 355
 356     if (overwrite)
 357     {
 358         mesh.setInstance(oldInstance);
 359     }
 360     Pout<< "Writing mesh to time " << runTime.timeName() << endl;
 361     mesh.write();
 362
 363     // Dump duplicated points (if any)
 364     if (split)
 365     {
 366         const labelList& pointMap = map().pointMap();
 367
 368         labelList nDupPerPoint(map().nOldPoints(), 0);
 369
 370         pointSet dupPoints(mesh, "duplicatedPoints", 100);
 371
 372         forAll(pointMap, pointI)
 373         {
 374             label oldPointI = pointMap[pointI];
 375
 376             nDupPerPoint[oldPointI]++;
 377
 378             if (nDupPerPoint[oldPointI] > 1)
 379             {
 380                 dupPoints.insert(map().reversePointMap()[oldPointI]);
 381                 dupPoints.insert(pointI);
 382             }
 383         }
 384
 385         Pout<< "Writing " << dupPoints.size()
 386             << " duplicated points to pointSet "
 387             << dupPoints.objectPath() << nl << endl;
 388
 389         dupPoints.write();
 390     }
 391
 392     Info<< "End\n" << endl;
 393
 394     return 0;
 395 }
 396
 397
 398 // ************************************************************************* //