applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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     couple
  44     {
  45         type                    slidingInterface;
  46         masterFaceZoneName      MSMasterZone
  47         slaveFaceZoneName       MSSlaveZone
  48         cutPointZoneName        MSCutPointZone
  49         cutFaceZoneName         MSCutFaceZone
  50         masterPatchName         M;
  51         slavePatchName          S;
  52         typeOfMatch             partial or integral
  53     }
  54
  55
  56 \*---------------------------------------------------------------------------*/
  57
  58 #include "fvCFD.H"
  59 #include "polyTopoChanger.H"
  60 #include "mapPolyMesh.H"
  61 #include "ListOps.H"
  62 #include "slidingInterface.H"
  63 #include "perfectInterface.H"
  64 #include "IOobjectList.H"
  65 #include "ReadFields.H"
  66
  67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  68
  69 // Checks whether patch present
  70 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
  71 {
  72     label patchI = bMesh.findPatchID(name);
  73
  74     if (patchI == -1)
  75     {
  76         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
  77             << "Cannot find patch " << name << endl
  78             << "It should be present and of non-zero size" << endl
  79             << "Valid patches are " << bMesh.names()
  80             << exit(FatalError);
  81     }
  82
  83     if (bMesh[patchI].empty())
  84     {
  85         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
  86             << "Patch " << name << " is present but zero size"
  87             << exit(FatalError);
  88     }
  89 }
  90
  91
  92 // Main program:
  93
  94 int main(int argc, char *argv[])
  95 {
  96     Foam::argList::noParallel();
  97 #   include "addRegionOption.H"
  98     Foam::argList::validArgs.append("masterPatch");
  99     Foam::argList::validArgs.append("slavePatch");
 100
 101     Foam::argList::validOptions.insert("partial", "");
 102     Foam::argList::validOptions.insert("perfect", "");
 103     Foam::argList::validOptions.insert("clearUnusedFaces", "");
 104     Foam::argList::validOptions.insert("overwrite", "");
 105
 106 #   include "setRootCase.H"
 107 #   include "createTime.H"
 108     runTime.functionObjects().off();
 109 #   include "createNamedMesh.H"
 110     const word oldInstance = mesh.pointsInstance();
 111
 112
 113     word masterPatchName(args.additionalArgs()[0]);
 114     word slavePatchName(args.additionalArgs()[1]);
 115
 116     bool partialCover = args.optionFound("partial");
 117     bool perfectCover = args.optionFound("perfect");
 118     bool clearUnusedFaces = args.optionFound("clearUnusedFaces");
 119     bool overwrite = args.optionFound("overwrite");
 120
 121     if (partialCover && perfectCover)
 122     {
 123         FatalErrorIn(args.executable())
 124             << "Cannot both supply partial and perfect." << endl
 125             << "Use perfect match option if the patches perfectly align"
 126             << " (both vertex positions and face centres)" << endl
 127             << exit(FatalError);
 128     }
 129
 130
 131     const word mergePatchName(masterPatchName + slavePatchName);
 132     const word cutZoneName(mergePatchName + "CutFaceZone");
 133
 134     slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
 135
 136     if (partialCover)
 137     {
 138         Info<< "Coupling partially overlapping patches "
 139             << masterPatchName << " and " << slavePatchName << nl
 140             << "Resulting internal faces will be in faceZone " << cutZoneName
 141             << nl
 142             << "Any uncovered faces will remain in their patch"
 143             << endl;
 144
 145         tom = slidingInterface::PARTIAL;
 146     }
 147     else if (perfectCover)
 148     {
 149         Info<< "Coupling perfectly aligned patches "
 150             << masterPatchName << " and " << slavePatchName << nl
 151             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 152             << nl << nl
 153             << "Note: both patches need to align perfectly." << nl
 154             << "Both the vertex"
 155             << " positions and the face centres need to align to within" << nl
 156             << "a tolerance given by the minimum edge length on the patch"
 157             << endl;
 158     }
 159     else
 160     {
 161         Info<< "Coupling patches " << masterPatchName << " and "
 162             << slavePatchName << nl
 163             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 164             << nl << nl
 165             << "Note: the overall area covered by both patches should be"
 166             << " identical (\"integral\" interface)." << endl
 167             << "If this is not the case use the -partial option" << nl << endl;
 168     }
 169
 170     // Check for non-empty master and slave patches
 171     checkPatch(mesh.boundaryMesh(), masterPatchName);
 172     checkPatch(mesh.boundaryMesh(), slavePatchName);
 173
 174     // Create and add face zones and mesh modifiers
 175
 176     // Master patch
 177     const polyPatch& masterPatch =
 178         mesh.boundaryMesh()
 179         [
 180             mesh.boundaryMesh().findPatchID(masterPatchName)
 181         ];
 182
 183     // Make list of masterPatch faces
 184     labelList isf(masterPatch.size());
 185
 186     forAll (isf, i)
 187     {
 188         isf[i] = masterPatch.start() + i;
 189     }
 190
 191     polyTopoChanger stitcher(mesh);
 192     stitcher.setSize(1);
 193
 194     DynamicList<pointZone*> pz;
 195     DynamicList<faceZone*> fz;
 196     DynamicList<cellZone*> cz;
 197
 198     if (perfectCover)
 199     {
 200         // Add empty zone for resulting internal faces
 201         fz.append
 202         (
 203             new faceZone
 204             (
 205                 cutZoneName,
 206                 isf,
 207                 boolList(masterPatch.size(), false),
 208                 0,
 209                 mesh.faceZones()
 210             )
 211         );
 212
 213         // Note: make sure to add the zones BEFORE constructing polyMeshModifier
 214         // (since looks up various zones at construction time)
 215         Info << "Adding point and face zones" << endl;
 216         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 217
 218         // Add the perfect interface mesh modifier
 219         stitcher.set
 220         (
 221             0,
 222             new perfectInterface
 223             (
 224                 "couple",
 225                 0,
 226                 stitcher,
 227                 cutZoneName,
 228                 masterPatchName,
 229                 slavePatchName
 230             )
 231         );
 232     }
 233     else
 234     {
 235         pz.append
 236         (
 237             new pointZone
 238             (
 239                 mergePatchName + "CutPointZone",
 240                 labelList(0),
 241                 0,
 242                 mesh.pointZones()
 243             )
 244         );
 245
 246         fz.append
 247         (
 248             new faceZone
 249             (
 250                 mergePatchName + "MasterZone",
 251                 isf,
 252                 boolList(masterPatch.size(), false),
 253                 0,
 254                 mesh.faceZones()
 255             )
 256         );
 257
 258         // Slave patch
 259         const polyPatch& slavePatch =
 260             mesh.boundaryMesh()
 261             [
 262                 mesh.boundaryMesh().findPatchID(slavePatchName)
 263             ];
 264
 265         labelList osf(slavePatch.size());
 266
 267         forAll (osf, i)
 268         {
 269             osf[i] = slavePatch.start() + i;
 270         }
 271
 272         fz.append
 273         (
 274             new faceZone
 275             (
 276                 mergePatchName + "SlaveZone",
 277                 osf,
 278                 boolList(slavePatch.size(), false),
 279                 1,
 280                 mesh.faceZones()
 281             )
 282         );
 283
 284         // Add empty zone for cut faces
 285         fz.append
 286         (
 287             new faceZone
 288             (
 289                 cutZoneName,
 290                 labelList(0),
 291                 boolList(0, false),
 292                 2,
 293                 mesh.faceZones()
 294             )
 295         );
 296
 297
 298         // Note: make sure to add the zones BEFORE constructing
 299         // polyMeshModifier (since looks up various zones at construction time)
 300         Info << "Adding point and face zones" << endl;
 301         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 302
 303         // Add the sliding interface mesh modifier
 304         stitcher.set
 305         (
 306             0,
 307             new slidingInterface
 308             (
 309                 "couple",
 310                 0,
 311                 stitcher,
 312                 mergePatchName + "MasterZone",
 313                 mergePatchName + "SlaveZone",
 314                 mergePatchName + "CutPointZone",
 315                 cutZoneName,
 316                 masterPatchName,
 317                 slavePatchName,
 318                 tom,                  // integral or partial
 319                 false,                // Attach-detach action
 320                 intersection::VISIBLE
 321             )
 322         );
 323     }
 324
 325
 326     // Search for list of objects for this time
 327     IOobjectList objects(mesh, runTime.timeName());
 328
 329     // Read all current fvFields so they will get mapped
 330     Info<< "Reading all current volfields" << endl;
 331     PtrList<volScalarField> volScalarFields;
 332     ReadFields(mesh, objects, volScalarFields);
 333
 334     PtrList<volVectorField> volVectorFields;
 335     ReadFields(mesh, objects, volVectorFields);
 336
 337     PtrList<volSphericalTensorField> volSphericalTensorFields;
 338     ReadFields(mesh, objects, volSphericalTensorFields);
 339
 340     PtrList<volSymmTensorField> volSymmTensorFields;
 341     ReadFields(mesh, objects, volSymmTensorFields);
 342
 343     PtrList<volTensorField> volTensorFields;
 344     ReadFields(mesh, objects, volTensorFields);
 345
 346     //- uncomment if you want to interpolate surface fields (usually bad idea)
 347     //Info<< "Reading all current surfaceFields" << endl;
 348     //PtrList<surfaceScalarField> surfaceScalarFields;
 349     //ReadFields(mesh, objects, surfaceScalarFields);
 350     //
 351     //PtrList<surfaceVectorField> surfaceVectorFields;
 352     //ReadFields(mesh, objects, surfaceVectorFields);
 353     //
 354     //PtrList<surfaceTensorField> surfaceTensorFields;
 355     //ReadFields(mesh, objects, surfaceTensorFields);
 356
 357     if (!overwrite)
 358     {
 359         runTime++;
 360     }
 361
 362     // Execute all polyMeshModifiers
 363     autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh();
 364
 365     mesh.movePoints(morphMap->preMotionPoints());
 366
 367     if (clearUnusedFaces)
 368     {
 369         // Clear unused points and faces by manually resetting the list"
 370         Info << "Clear unused points and faces" << nl << endl;
 371
 372         pointField& p = const_cast<pointField&>(mesh.allPoints());
 373         p.setSize(mesh.nPoints());
 374
 375         faceList& f = const_cast<faceList&>(mesh.allFaces());
 376         f.setSize(mesh.nFaces());
 377
 378         Xfer<pointField> pXfer(p);
 379         Xfer<faceList> fXfer(f);
 380         Xfer<labelList> ownXfer(mesh.faceOwner());
 381         Xfer<labelList> neiXfer(mesh.faceNeighbour());
 382
 383         label nPatches = mesh.boundaryMesh().size();
 384         labelList patchSizes(nPatches, 0);
 385         labelList patchStarts(nPatches, -1);
 386         forAll(patchSizes, patchI)
 387         {
 388             patchSizes[patchI] = mesh.boundaryMesh()[patchI].size();
 389             patchStarts[patchI] = mesh.boundaryMesh()[patchI].start();
 390         }
 391
 392         mesh.removeZones();
 393
 394         mesh.resetPrimitives
 395         (
 396             pXfer,
 397             fXfer,
 398             ownXfer,
 399             neiXfer,
 400             patchSizes,
 401             patchStarts
 402         );
 403     }
 404
 405
 406     // Write mesh
 407     if (overwrite)
 408     {
 409         mesh.setInstance(oldInstance);
 410         stitcher.instance() = oldInstance;
 411     }
 412     Info << nl << "Writing polyMesh to time " << runTime.timeName() << endl;
 413
 414     IOstream::defaultPrecision(10);
 415
 416     // Bypass runTime write (since only writes at outputTime)
 417     if
 418     (
 419        !runTime.objectRegistry::writeObject
 420         (
 421             runTime.writeFormat(),
 422             IOstream::currentVersion,
 423             runTime.writeCompression()
 424         )
 425     )
 426     {
 427         FatalErrorIn(args.executable())
 428             << "Failed writing polyMesh."
 429             << exit(FatalError);
 430     }
 431
 432     mesh.faceZones().write();
 433     mesh.pointZones().write();
 434     mesh.cellZones().write();
 435
 436     // Write fields
 437     runTime.write();
 438
 439     Info<< nl << "end" << endl;
 440
 441     return 0;
 442 }
 443
 444
 445 // ************************************************************************* //