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) 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     'Stitches' a mesh.
  27
  28     Takes a mesh and two patches and merges the faces on the two patches
  29     (if geometrically possible) so the faces become internal.
  30
  31     Can do
  32     - 'perfect' match: faces and points on patches align exactly. Order might
  33     be different though.
  34     - 'integral' match: where the surfaces on both patches exactly
  35     match but the individual faces not
  36     - 'partial' match: where the non-overlapping part of the surface remains
  37     in the respective patch.
  38
  39     Note : Is just a front-end to perfectInterface/slidingInterface.
  40
  41     Comparable to running a meshModifier of the form
  42     (if masterPatch is called "M" and slavePatch "S"):
  43
  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
  56
  57 \*---------------------------------------------------------------------------*/
  58
  59 #include "fvCFD.H"
  60 #include "polyTopoChanger.H"
  61 #include "mapPolyMesh.H"
  62 #include "ListOps.H"
  63 #include "slidingInterface.H"
  64 #include "perfectInterface.H"
  65 #include "IOobjectList.H"
  66 #include "ReadFields.H"
  67
  68
  69 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  70
  71 // Checks whether patch present
  72 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
  73 {
  74     label patchI = bMesh.findPatchID(name);
  75
  76     if (patchI == -1)
  77     {
  78         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
  79             << "Cannot find patch " << name << endl
  80             << "It should be present and of non-zero size" << endl
  81             << "Valid patches are " << bMesh.names()
  82             << exit(FatalError);
  83     }
  84
  85     if (bMesh[patchI].size() == 0)
  86     {
  87         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
  88             << "Patch " << name << " is present but zero size"
  89             << exit(FatalError);
  90     }
  91 }
  92
  93
  94 // Main program:
  95
  96 int main(int argc, char *argv[])
  97 {
  98     Foam::argList::noParallel();
  99
 100     Foam::argList::validArgs.append("masterPatch");
 101     Foam::argList::validArgs.append("slavePatch");
 102
 103     Foam::argList::validOptions.insert("partial", "");
 104     Foam::argList::validOptions.insert("perfect", "");
 105
 106     Foam::argList::validOptions.insert("overwrite", "");
 107
 108 #   include "setRootCase.H"
 109 #   include "createTime.H"
 110     runTime.functionObjects().off();
 111 #   include "createMesh.H"
 112     const word oldInstance = mesh.pointsInstance();
 113
 114
 115     word masterPatchName(args.additionalArgs()[0]);
 116     word slavePatchName(args.additionalArgs()[1]);
 117
 118     bool partialCover = args.options().found("partial");
 119     bool perfectCover = args.options().found("perfect");
 120     bool overwrite = args.options().found("overwrite");
 121
 122     if (partialCover && perfectCover)
 123     {
 124         FatalErrorIn(args.executable())
 125             << "Cannot both supply partial and perfect." << endl
 126             << "Use perfect match option if the patches perfectly align"
 127             << " (both vertex positions and face centres)" << endl
 128             << exit(FatalError);
 129     }
 130
 131
 132     const word mergePatchName(masterPatchName + slavePatchName);
 133     const word cutZoneName(mergePatchName + "CutFaceZone");
 134
 135     slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
 136
 137     if (partialCover)
 138     {
 139         Info<< "Coupling partially overlapping patches "
 140             << masterPatchName << " and " << slavePatchName << nl
 141             << "Resulting internal faces will be in faceZone " << cutZoneName
 142             << nl
 143             << "Any uncovered faces will remain in their patch"
 144             << endl;
 145
 146         tom = slidingInterface::PARTIAL;
 147     }
 148     else if (perfectCover)
 149     {
 150         Info<< "Coupling perfectly aligned patches "
 151             << masterPatchName << " and " << slavePatchName << nl
 152             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 153             << nl << nl
 154             << "Note: both patches need to align perfectly." << nl
 155             << "Both the vertex"
 156             << " positions and the face centres need to align to within" << nl
 157             << "a tolerance given by the minimum edge length on the patch"
 158             << endl;
 159     }
 160     else
 161     {
 162         Info<< "Coupling patches " << masterPatchName << " and "
 163             << slavePatchName << nl
 164             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 165             << nl << nl
 166             << "Note: the overall area covered by both patches should be"
 167             << " identical (\"integral\" interface)." << endl
 168             << "If this is not the case use the -partial option" << nl << endl;
 169     }
 170
 171     // Check for non-empty master and slave patches
 172     checkPatch(mesh.boundaryMesh(), masterPatchName);
 173     checkPatch(mesh.boundaryMesh(), slavePatchName);
 174
 175     // Create and add face zones and mesh modifiers
 176
 177     // Master patch
 178     const polyPatch& masterPatch =
 179         mesh.boundaryMesh()
 180         [
 181             mesh.boundaryMesh().findPatchID(masterPatchName)
 182         ];
 183
 184     // Make list of masterPatch faces
 185     labelList isf(masterPatch.size());
 186
 187     forAll (isf, i)
 188     {
 189         isf[i] = masterPatch.start() + i;
 190     }
 191
 192     polyTopoChanger stitcher(mesh);
 193     stitcher.setSize(1);
 194
 195     DynamicList<pointZone*> pz;
 196     DynamicList<faceZone*> fz;
 197     DynamicList<cellZone*> cz;
 198
 199     if (perfectCover)
 200     {
 201         // Add empty zone for resulting internal faces
 202         fz.append
 203         (
 204             new faceZone
 205             (
 206                 cutZoneName,
 207                 isf,
 208                 boolList(masterPatch.size(), false),
 209                 0,
 210                 mesh.faceZones()
 211             )
 212         );
 213
 214         // Note: make sure to add the zones BEFORE constructing polyMeshModifier
 215         // (since looks up various zones at construction time)
 216         Info << "Adding point and face zones" << endl;
 217         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 218
 219         // Add the perfect interface mesh modifier
 220         stitcher.set
 221         (
 222             0,
 223             new perfectInterface
 224             (
 225                 "couple",
 226                 0,
 227                 stitcher,
 228                 cutZoneName,
 229                 masterPatchName,
 230                 slavePatchName
 231             )
 232         );
 233     }
 234     else
 235     {
 236         pz.append
 237         (
 238             new pointZone
 239             (
 240                 mergePatchName + "CutPointZone",
 241                 labelList(0),
 242                 0,
 243                 mesh.pointZones()
 244             )
 245         );
 246
 247         fz.append
 248         (
 249             new faceZone
 250             (
 251                 mergePatchName + "MasterZone",
 252                 isf,
 253                 boolList(masterPatch.size(), false),
 254                 0,
 255                 mesh.faceZones()
 256             )
 257         );
 258
 259         // Slave patch
 260         const polyPatch& slavePatch =
 261             mesh.boundaryMesh()
 262             [
 263                 mesh.boundaryMesh().findPatchID(slavePatchName)
 264             ];
 265
 266         labelList osf(slavePatch.size());
 267
 268         forAll (osf, i)
 269         {
 270             osf[i] = slavePatch.start() + i;
 271         }
 272
 273         fz.append
 274         (
 275             new faceZone
 276             (
 277                 mergePatchName + "SlaveZone",
 278                 osf,
 279                 boolList(slavePatch.size(), false),
 280                 1,
 281                 mesh.faceZones()
 282             )
 283         );
 284
 285         // Add empty zone for cut faces
 286         fz.append
 287         (
 288             new faceZone
 289             (
 290                 cutZoneName,
 291                 labelList(0),
 292                 boolList(0, false),
 293                 2,
 294                 mesh.faceZones()
 295             )
 296         );
 297
 298
 299         // Note: make sure to add the zones BEFORE constructing polyMeshModifier
 300         // (since looks up various zones at construction time)
 301         Info << "Adding point and face zones" << endl;
 302         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 303
 304         // Add the sliding interface mesh modifier
 305         stitcher.set
 306         (
 307             0,
 308             new slidingInterface
 309             (
 310                 "couple",
 311                 0,
 312                 stitcher,
 313                 mergePatchName + "MasterZone",
 314                 mergePatchName + "SlaveZone",
 315                 mergePatchName + "CutPointZone",
 316                 cutZoneName,
 317                 masterPatchName,
 318                 slavePatchName,
 319                 tom,                    // integral or partial
 320                 true                    // couple/decouple mode
 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(true);
 364
 365     mesh.movePoints(morphMap->preMotionPoints());
 366
 367     // Write mesh
 368     if (overwrite)
 369     {
 370         mesh.setInstance(oldInstance);
 371     }
 372     Info << nl << "Writing polyMesh to time " << runTime.timeName() << endl;
 373
 374     IOstream::defaultPrecision(10);
 375
 376     // Bypass runTime write (since only writes at outputTime)
 377     if
 378     (
 379        !runTime.objectRegistry::writeObject
 380         (
 381             runTime.writeFormat(),
 382             IOstream::currentVersion,
 383             runTime.writeCompression()
 384         )
 385     )
 386     {
 387         FatalErrorIn(args.executable())
 388             << "Failed writing polyMesh."
 389             << exit(FatalError);
 390     }
 391
 392     mesh.faceZones().write();
 393     mesh.pointZones().write();
 394     mesh.cellZones().write();
 395
 396     // Write fields
 397     runTime.write();
 398
 399     Info<< nl << "end" << endl;
 400
 401     return 0;
 402 }
 403
 404
 405 // ************************************************************************* //