applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.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     'Stitches' a mesh.
  26
  27     Takes a mesh and two patches and merges the faces on the two patches
  28     (if geometrically possible) so the faces become internal.
  29
  30     Can do
  31     - 'perfect' match: faces and points on patches align exactly. Order might
  32     be different though.
  33     - 'integral' match: where the surfaces on both patches exactly
  34     match but the individual faces not
  35     - 'partial' match: where the non-overlapping part of the surface remains
  36     in the respective patch.
  37
  38     Note : Is just a front-end to perfectInterface/slidingInterface.
  39
  40     Comparable to running a meshModifier of the form
  41     (if masterPatch is called "M" and slavePatch "S"):
  42
  43     \verbatim
  44     couple
  45     {
  46         type                    slidingInterface;
  47         masterFaceZoneName      MSMasterZone
  48         slaveFaceZoneName       MSSlaveZone
  49         cutPointZoneName        MSCutPointZone
  50         cutFaceZoneName         MSCutFaceZone
  51         masterPatchName         M;
  52         slavePatchName          S;
  53         typeOfMatch             partial or integral
  54     }
  55     \endverbatim
  56
  57
  58 \*---------------------------------------------------------------------------*/
  59
  60 #include "fvCFD.H"
  61 #include "polyTopoChanger.H"
  62 #include "mapPolyMesh.H"
  63 #include "ListOps.H"
  64 #include "slidingInterface.H"
  65 #include "perfectInterface.H"
  66 #include "IOobjectList.H"
  67 #include "ReadFields.H"
  68
  69
  70 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  71
  72 label addPointZone(const polyMesh& mesh, const word& name)
  73 {
  74     label zoneID = mesh.pointZones().findZoneID(name);
  75
  76     if (zoneID != -1)
  77     {
  78         Info<< "Reusing existing pointZone "
  79             << mesh.pointZones()[zoneID].name()
  80             << " at index " << zoneID << endl;
  81     }
  82     else
  83     {
  84         pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh).pointZones();
  85         zoneID = pointZones.size();
  86         Info<< "Adding pointZone " << name << " at index " << zoneID << endl;
  87
  88         pointZones.setSize(zoneID+1);
  89         pointZones.set
  90         (
  91             zoneID,
  92             new pointZone
  93             (
  94                 name,
  95                 labelList(0),
  96                 zoneID,
  97                 pointZones
  98             )
  99         );
 100     }
 101     return zoneID;
 102 }
 103
 104
 105 label addFaceZone(const polyMesh& mesh, const word& name)
 106 {
 107     label zoneID = mesh.faceZones().findZoneID(name);
 108
 109     if (zoneID != -1)
 110     {
 111         Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
 112             << " at index " << zoneID << endl;
 113     }
 114     else
 115     {
 116         faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh).faceZones();
 117         zoneID = faceZones.size();
 118         Info<< "Adding faceZone " << name << " at index " << zoneID << endl;
 119
 120         faceZones.setSize(zoneID+1);
 121         faceZones.set
 122         (
 123             zoneID,
 124             new faceZone
 125             (
 126                 name,
 127                 labelList(0),
 128                 boolList(),
 129                 zoneID,
 130                 faceZones
 131             )
 132         );
 133     }
 134     return zoneID;
 135 }
 136
 137
 138 label addCellZone(const polyMesh& mesh, const word& name)
 139 {
 140     label zoneID = mesh.cellZones().findZoneID(name);
 141
 142     if (zoneID != -1)
 143     {
 144         Info<< "Reusing existing cellZone " << mesh.cellZones()[zoneID].name()
 145             << " at index " << zoneID << endl;
 146     }
 147     else
 148     {
 149         cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh).cellZones();
 150         zoneID = cellZones.size();
 151         Info<< "Adding cellZone " << name << " at index " << zoneID << endl;
 152
 153         cellZones.setSize(zoneID+1);
 154         cellZones.set
 155         (
 156             zoneID,
 157             new cellZone
 158             (
 159                 name,
 160                 labelList(0),
 161                 zoneID,
 162                 cellZones
 163             )
 164         );
 165     }
 166     return zoneID;
 167 }
 168
 169
 170 // Checks whether patch present
 171 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
 172 {
 173     const label patchI = bMesh.findPatchID(name);
 174
 175     if (patchI == -1)
 176     {
 177         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
 178             << "Cannot find patch " << name << endl
 179             << "It should be present and of non-zero size" << endl
 180             << "Valid patches are " << bMesh.names()
 181             << exit(FatalError);
 182     }
 183
 184     if (bMesh[patchI].empty())
 185     {
 186         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
 187             << "Patch " << name << " is present but zero size"
 188             << exit(FatalError);
 189     }
 190 }
 191
 192
 193 // Main program:
 194
 195 int main(int argc, char *argv[])
 196 {
 197     argList::addNote
 198     (
 199         "merge the faces on the specified patches (if geometrically possible)\n"
 200         "so the faces become internal"
 201     );
 202
 203     argList::noParallel();
 204     #include "addOverwriteOption.H"
 205     #include "addRegionOption.H"
 206
 207     argList::validArgs.append("masterPatch");
 208     argList::validArgs.append("slavePatch");
 209
 210     argList::addBoolOption
 211     (
 212         "partial",
 213         "couple partially overlapping patches"
 214     );
 215     argList::addBoolOption
 216     (
 217         "perfect",
 218         "couple perfectly aligned patches"
 219     );
 220     argList::addOption
 221     (
 222         "toleranceDict",
 223         "file",
 224         "dictionary file with tolerances"
 225     );
 226
 227     #include "setRootCase.H"
 228     #include "createTime.H"
 229     runTime.functionObjects().off();
 230     #include "createNamedMesh.H"
 231
 232     const word oldInstance = mesh.pointsInstance();
 233
 234     const word masterPatchName = args[1];
 235     const word slavePatchName  = args[2];
 236
 237     const bool partialCover = args.optionFound("partial");
 238     const bool perfectCover = args.optionFound("perfect");
 239     const bool overwrite    = args.optionFound("overwrite");
 240
 241     if (partialCover && perfectCover)
 242     {
 243         FatalErrorIn(args.executable())
 244             << "Cannot supply both partial and perfect." << endl
 245             << "Use perfect match option if the patches perfectly align"
 246             << " (both vertex positions and face centres)" << endl
 247             << exit(FatalError);
 248     }
 249
 250
 251     const word mergePatchName(masterPatchName + slavePatchName);
 252     const word cutZoneName(mergePatchName + "CutFaceZone");
 253
 254     slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
 255
 256     if (partialCover)
 257     {
 258         Info<< "Coupling partially overlapping patches "
 259             << masterPatchName << " and " << slavePatchName << nl
 260             << "Resulting internal faces will be in faceZone " << cutZoneName
 261             << nl
 262             << "Any uncovered faces will remain in their patch"
 263             << endl;
 264
 265         tom = slidingInterface::PARTIAL;
 266     }
 267     else if (perfectCover)
 268     {
 269         Info<< "Coupling perfectly aligned patches "
 270             << masterPatchName << " and " << slavePatchName << nl
 271             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 272             << nl << nl
 273             << "Note: both patches need to align perfectly." << nl
 274             << "Both the vertex"
 275             << " positions and the face centres need to align to within" << nl
 276             << "a tolerance given by the minimum edge length on the patch"
 277             << endl;
 278     }
 279     else
 280     {
 281         Info<< "Coupling patches " << masterPatchName << " and "
 282             << slavePatchName << nl
 283             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 284             << nl << nl
 285             << "Note: the overall area covered by both patches should be"
 286             << " identical (\"integral\" interface)." << endl
 287             << "If this is not the case use the -partial option" << nl << endl;
 288     }
 289
 290     // set up the tolerances for the sliding mesh
 291     dictionary slidingTolerances;
 292     if (args.options().found("toleranceDict"))
 293     {
 294         IOdictionary toleranceFile
 295         (
 296             IOobject
 297             (
 298                 args.options()["toleranceDict"],
 299                 runTime.constant(),
 300                 mesh,
 301                 IOobject::MUST_READ_IF_MODIFIED,
 302                 IOobject::NO_WRITE
 303             )
 304         );
 305         slidingTolerances += toleranceFile;
 306     }
 307
 308     // Check for non-empty master and slave patches
 309     checkPatch(mesh.boundaryMesh(), masterPatchName);
 310     checkPatch(mesh.boundaryMesh(), slavePatchName);
 311
 312     // Create and add face zones and mesh modifiers
 313
 314     // Master patch
 315     const polyPatch& masterPatch = mesh.boundaryMesh()[masterPatchName];
 316
 317     // Make list of masterPatch faces
 318     labelList isf(masterPatch.size());
 319
 320     forAll(isf, i)
 321     {
 322         isf[i] = masterPatch.start() + i;
 323     }
 324
 325     polyTopoChanger stitcher(mesh);
 326     stitcher.setSize(1);
 327
 328     mesh.pointZones().clearAddressing();
 329     mesh.faceZones().clearAddressing();
 330     mesh.cellZones().clearAddressing();
 331
 332     if (perfectCover)
 333     {
 334         // Add empty zone for resulting internal faces
 335         label cutZoneID = addFaceZone(mesh, cutZoneName);
 336
 337         mesh.faceZones()[cutZoneID].resetAddressing
 338         (
 339             isf,
 340             boolList(masterPatch.size(), false)
 341         );
 342
 343         // Add the perfect interface mesh modifier
 344         stitcher.set
 345         (
 346             0,
 347             new perfectInterface
 348             (
 349                 "couple",
 350                 0,
 351                 stitcher,
 352                 cutZoneName,
 353                 masterPatchName,
 354                 slavePatchName
 355             )
 356         );
 357     }
 358     else
 359     {
 360         label pointZoneID = addPointZone(mesh, mergePatchName + "CutPointZone");
 361         mesh.pointZones()[pointZoneID] = labelList(0);
 362
 363         label masterZoneID = addFaceZone(mesh, mergePatchName + "MasterZone");
 364
 365         mesh.faceZones()[masterZoneID].resetAddressing
 366         (
 367             isf,
 368             boolList(masterPatch.size(), false)
 369         );
 370
 371         // Slave patch
 372         const polyPatch& slavePatch = mesh.boundaryMesh()[slavePatchName];
 373
 374         labelList osf(slavePatch.size());
 375
 376         forAll(osf, i)
 377         {
 378             osf[i] = slavePatch.start() + i;
 379         }
 380
 381         label slaveZoneID = addFaceZone(mesh, mergePatchName + "SlaveZone");
 382         mesh.faceZones()[slaveZoneID].resetAddressing
 383         (
 384             osf,
 385             boolList(slavePatch.size(), false)
 386         );
 387
 388         // Add empty zone for cut faces
 389         label cutZoneID = addFaceZone(mesh, cutZoneName);
 390         mesh.faceZones()[cutZoneID].resetAddressing
 391         (
 392             labelList(0),
 393             boolList(0, false)
 394         );
 395
 396
 397         // Add the sliding interface mesh modifier
 398         stitcher.set
 399         (
 400             0,
 401             new slidingInterface
 402             (
 403                 "couple",
 404                 0,
 405                 stitcher,
 406                 mergePatchName + "MasterZone",
 407                 mergePatchName + "SlaveZone",
 408                 mergePatchName + "CutPointZone",
 409                 cutZoneName,
 410                 masterPatchName,
 411                 slavePatchName,
 412                 tom,                    // integral or partial
 413                 true                    // couple/decouple mode
 414             )
 415         );
 416         static_cast<slidingInterface&>(stitcher[0]).setTolerances
 417         (
 418             slidingTolerances,
 419             true
 420         );
 421     }
 422
 423
 424     // Search for list of objects for this time
 425     IOobjectList objects(mesh, runTime.timeName());
 426
 427     // Read all current fvFields so they will get mapped
 428     Info<< "Reading all current volfields" << endl;
 429     PtrList<volScalarField> volScalarFields;
 430     ReadFields(mesh, objects, volScalarFields);
 431
 432     PtrList<volVectorField> volVectorFields;
 433     ReadFields(mesh, objects, volVectorFields);
 434
 435     PtrList<volSphericalTensorField> volSphericalTensorFields;
 436     ReadFields(mesh, objects, volSphericalTensorFields);
 437
 438     PtrList<volSymmTensorField> volSymmTensorFields;
 439     ReadFields(mesh, objects, volSymmTensorFields);
 440
 441     PtrList<volTensorField> volTensorFields;
 442     ReadFields(mesh, objects, volTensorFields);
 443
 444     //- uncomment if you want to interpolate surface fields (usually bad idea)
 445     //Info<< "Reading all current surfaceFields" << endl;
 446     //PtrList<surfaceScalarField> surfaceScalarFields;
 447     //ReadFields(mesh, objects, surfaceScalarFields);
 448     //
 449     //PtrList<surfaceVectorField> surfaceVectorFields;
 450     //ReadFields(mesh, objects, surfaceVectorFields);
 451     //
 452     //PtrList<surfaceTensorField> surfaceTensorFields;
 453     //ReadFields(mesh, objects, surfaceTensorFields);
 454
 455     if (!overwrite)
 456     {
 457         runTime++;
 458     }
 459
 460     // Execute all polyMeshModifiers
 461     autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
 462
 463     mesh.movePoints(morphMap->preMotionPoints());
 464
 465     // Write mesh
 466     if (overwrite)
 467     {
 468         mesh.setInstance(oldInstance);
 469         stitcher.instance() = oldInstance;
 470     }
 471     Info<< nl << "Writing polyMesh to time " << runTime.timeName() << endl;
 472
 473     IOstream::defaultPrecision(10);
 474
 475     // Bypass runTime write (since only writes at outputTime)
 476     if
 477     (
 478        !runTime.objectRegistry::writeObject
 479         (
 480             runTime.writeFormat(),
 481             IOstream::currentVersion,
 482             runTime.writeCompression()
 483         )
 484     )
 485     {
 486         FatalErrorIn(args.executable())
 487             << "Failed writing polyMesh."
 488             << exit(FatalError);
 489     }
 490
 491     mesh.faceZones().write();
 492     mesh.pointZones().write();
 493     mesh.cellZones().write();
 494
 495     // Write fields
 496     runTime.write();
 497
 498     Info<< nl << "end" << endl;
 499
 500     return 0;
 501 }
 502
 503
 504 // ************************************************************************* //