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-2008 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 "IndirectList.H"
  64 #include "slidingInterface.H"
  65 #include "perfectInterface.H"
  66 #include "IOobjectList.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 // Read field
  95 template<class GeoField>
  96 void readFields
  97 (
  98     const fvMesh& mesh,
  99     const IOobjectList& objects,
 100     PtrList<GeoField>& fields
 101 )
 102 {
 103     // Search list of objects for volScalarFields
 104     IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
 105
 106     // Construct the vol scalar fields
 107     fields.setSize(fieldObjects.size());
 108
 109     label fieldi = 0;
 110     for
 111     (
 112         IOobjectList::iterator iter = fieldObjects.begin();
 113         iter != fieldObjects.end();
 114         ++iter
 115     )
 116     {
 117         fields.set(fieldi++, new GeoField(*iter(), mesh));
 118     }
 119 }
 120
 121
 122 // Main program:
 123
 124 int main(int argc, char *argv[])
 125 {
 126     Foam::argList::noParallel();
 127
 128     Foam::argList::validArgs.append("masterPatch");
 129     Foam::argList::validArgs.append("slavePatch");
 130
 131     Foam::argList::validOptions.insert("partial", "");
 132     Foam::argList::validOptions.insert("perfect", "");
 133
 134     Foam::argList::validOptions.insert("overwrite", "");
 135
 136 #   include "setRootCase.H"
 137 #   include "createTime.H"
 138 #   include "createMesh.H"
 139
 140
 141     word masterPatchName(args.additionalArgs()[0]);
 142     word slavePatchName(args.additionalArgs()[1]);
 143
 144     bool partialCover = args.options().found("partial");
 145     bool perfectCover = args.options().found("perfect");
 146     bool overwrite = args.options().found("overwrite");
 147
 148     if (partialCover && perfectCover)
 149     {
 150         FatalErrorIn(args.executable())
 151             << "Cannot both supply partial and perfect." << endl
 152             << "Use perfect match option if the patches perfectly align"
 153             << " (both vertex positions and face centres)" << endl
 154             << exit(FatalError);
 155     }
 156
 157
 158     const word mergePatchName(masterPatchName + slavePatchName);
 159     const word cutZoneName(mergePatchName + "CutFaceZone");
 160
 161     slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
 162
 163     if (partialCover)
 164     {
 165         Info<< "Coupling partially overlapping patches "
 166             << masterPatchName << " and " << slavePatchName << nl
 167             << "Resulting internal faces will be in faceZone " << cutZoneName
 168             << nl
 169             << "Any uncovered faces will remain in their patch"
 170             << endl;
 171
 172         tom = slidingInterface::PARTIAL;
 173     }
 174     else if (perfectCover)
 175     {
 176         Info<< "Coupling perfectly aligned patches "
 177             << masterPatchName << " and " << slavePatchName << nl
 178             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 179             << nl << nl
 180             << "Note: both patches need to align perfectly." << nl
 181             << "Both the vertex"
 182             << " positions and the face centres need to align to within" << nl
 183             << "a tolerance given by the minimum edge length on the patch"
 184             << endl;
 185     }
 186     else
 187     {
 188         Info<< "Coupling patches " << masterPatchName << " and "
 189             << slavePatchName << nl
 190             << "Resulting (internal) faces will be in faceZone " << cutZoneName
 191             << nl << nl
 192             << "Note: the overall area covered by both patches should be"
 193             << " identical (\"integral\" interface)." << endl
 194             << "If this is not the case use the -partial option" << nl << endl;
 195     }
 196
 197     // Check for non-empty master and slave patches
 198     checkPatch(mesh.boundaryMesh(), masterPatchName);
 199     checkPatch(mesh.boundaryMesh(), slavePatchName);
 200
 201     // Create and add face zones and mesh modifiers
 202
 203     // Master patch
 204     const polyPatch& masterPatch =
 205         mesh.boundaryMesh()
 206         [
 207             mesh.boundaryMesh().findPatchID(masterPatchName)
 208         ];
 209
 210     // Make list of masterPatch faces
 211     labelList isf(masterPatch.size());
 212
 213     forAll (isf, i)
 214     {
 215         isf[i] = masterPatch.start() + i;
 216     }
 217
 218     polyTopoChanger stitcher(mesh);
 219     stitcher.setSize(1);
 220
 221     DynamicList<pointZone*> pz;
 222     DynamicList<faceZone*> fz;
 223     DynamicList<cellZone*> cz;
 224
 225     if (perfectCover)
 226     {
 227         // Add empty zone for resulting internal faces
 228         fz.append
 229         (
 230             new faceZone
 231             (
 232                 cutZoneName,
 233                 isf,
 234                 boolList(masterPatch.size(), false),
 235                 0,
 236                 mesh.faceZones()
 237             )
 238         );
 239
 240         // Note: make sure to add the zones BEFORE constructing polyMeshModifier
 241         // (since looks up various zones at construction time)
 242         Info << "Adding point and face zones" << endl;
 243         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 244
 245         // Add the perfect interface mesh modifier
 246         stitcher.set
 247         (
 248             0,
 249             new perfectInterface
 250             (
 251                 "couple",
 252                 0,
 253                 stitcher,
 254                 cutZoneName,
 255                 masterPatchName,
 256                 slavePatchName
 257             )
 258         );
 259     }
 260     else
 261     {
 262         pz.append
 263         (
 264             new pointZone
 265             (
 266                 mergePatchName + "CutPointZone",
 267                 labelList(0),
 268                 0,
 269                 mesh.pointZones()
 270             )
 271         );
 272
 273         fz.append
 274         (
 275             new faceZone
 276             (
 277                 mergePatchName + "MasterZone",
 278                 isf,
 279                 boolList(masterPatch.size(), false),
 280                 0,
 281                 mesh.faceZones()
 282             )
 283         );
 284
 285         // Slave patch
 286         const polyPatch& slavePatch =
 287             mesh.boundaryMesh()
 288             [
 289                 mesh.boundaryMesh().findPatchID(slavePatchName)
 290             ];
 291
 292         labelList osf(slavePatch.size());
 293
 294         forAll (osf, i)
 295         {
 296             osf[i] = slavePatch.start() + i;
 297         }
 298
 299         fz.append
 300         (
 301             new faceZone
 302             (
 303                 mergePatchName + "SlaveZone",
 304                 osf,
 305                 boolList(slavePatch.size(), false),
 306                 1,
 307                 mesh.faceZones()
 308             )
 309         );
 310
 311         // Add empty zone for cut faces
 312         fz.append
 313         (
 314             new faceZone
 315             (
 316                 cutZoneName,
 317                 labelList(0),
 318                 boolList(0, false),
 319                 2,
 320                 mesh.faceZones()
 321             )
 322         );
 323
 324
 325         // Note: make sure to add the zones BEFORE constructing polyMeshModifier
 326         // (since looks up various zones at construction time)
 327         Info << "Adding point and face zones" << endl;
 328         mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
 329
 330         // Add the sliding interface mesh modifier
 331         stitcher.set
 332         (
 333             0,
 334             new slidingInterface
 335             (
 336                 "couple",
 337                 0,
 338                 stitcher,
 339                 mergePatchName + "MasterZone",
 340                 mergePatchName + "SlaveZone",
 341                 mergePatchName + "CutPointZone",
 342                 cutZoneName,
 343                 masterPatchName,
 344                 slavePatchName,
 345                 tom                   // integral or partial
 346             )
 347         );
 348     }
 349
 350
 351     // Search for list of objects for this time
 352     IOobjectList objects(mesh, runTime.timeName());
 353
 354     // Read all current fvFields so they will get mapped
 355     Info<< "Reading all current volfields" << endl;
 356     PtrList<volScalarField> volScalarFields;
 357     readFields(mesh, objects, volScalarFields);
 358
 359     PtrList<volVectorField> volVectorFields;
 360     readFields(mesh, objects, volVectorFields);
 361
 362     PtrList<volSphericalTensorField> volSphericalTensorFields;
 363     readFields(mesh, objects, volSphericalTensorFields);
 364
 365     PtrList<volSymmTensorField> volSymmTensorFields;
 366     readFields(mesh, objects, volSymmTensorFields);
 367
 368     PtrList<volTensorField> volTensorFields;
 369     readFields(mesh, objects, volTensorFields);
 370
 371     //- uncomment if you want to interpolate surface fields (usually bad idea)
 372     //Info<< "Reading all current surfaceFields" << endl;
 373     //PtrList<surfaceScalarField> surfaceScalarFields;
 374     //readFields(mesh, objects, surfaceScalarFields);
 375     //
 376     //PtrList<surfaceVectorField> surfaceVectorFields;
 377     //readFields(mesh, objects, surfaceVectorFields);
 378     //
 379     //PtrList<surfaceTensorField> surfaceTensorFields;
 380     //readFields(mesh, objects, surfaceTensorFields);
 381
 382     if (!overwrite)
 383     {
 384         runTime++;
 385     }
 386
 387     // Execute all polyMeshModifiers
 388     autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
 389
 390     mesh.movePoints(morphMap->preMotionPoints());
 391
 392     // Write mesh
 393     Info << nl << "Writing polyMesh to time " << runTime.timeName() << endl;
 394
 395     IOstream::defaultPrecision(10);
 396     if (!mesh.write())
 397     {
 398         FatalErrorIn(args.executable())
 399             << "Failed writing polyMesh."
 400             << exit(FatalError);
 401     }
 402
 403     mesh.faceZones().write();
 404     mesh.pointZones().write();
 405     mesh.cellZones().write();
 406
 407     // Write fields
 408     runTime.write();
 409
 410     Info<< nl << "end" << endl;
 411
 412     return 0;
 413 }
 414
 415
 416 // ************************************************************************* //