applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.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     Selects a section of mesh based on a cellSet.
  27
  28     The utility sub-sets the mesh to choose only a part of interest. Check
  29     the setSet/cellSet utilities to see how to select cells based on various.
  30
  31     The mesh will subset all points, faces and cells needed to make a sub-mesh
  32     but will not preserve attached boundary types.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "fvMeshSubset.H"
  37 #include "argList.H"
  38 #include "cellSet.H"
  39 #include "IOobjectList.H"
  40 #include "volFields.H"
  41
  42 using namespace Foam;
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46
  47 template<class Type>
  48 void subsetVolFields
  49 (
  50     const fvMeshSubset& subsetter,
  51     const wordList& fieldNames,
  52     PtrList<GeometricField<Type, fvPatchField, volMesh> >& subFields
  53 )
  54 {
  55     const fvMesh& baseMesh = subsetter.baseMesh();
  56
  57     forAll(fieldNames, i)
  58     {
  59         const word& fieldName = fieldNames[i];
  60
  61         Info<< "Subsetting field " << fieldName << endl;
  62
  63         GeometricField<Type, fvPatchField, volMesh> fld
  64         (
  65             IOobject
  66             (
  67                 fieldName,
  68                 baseMesh.time().timeName(),
  69                 baseMesh,
  70                 IOobject::MUST_READ,
  71                 IOobject::NO_WRITE
  72             ),
  73             baseMesh
  74         );
  75
  76         subFields.set(i, subsetter.interpolate(fld));
  77     }
  78 }
  79
  80
  81 template<class Type>
  82 void subsetSurfaceFields
  83 (
  84     const fvMeshSubset& subsetter,
  85     const wordList& fieldNames,
  86     PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& subFields
  87 )
  88 {
  89     const fvMesh& baseMesh = subsetter.baseMesh();
  90
  91     forAll(fieldNames, i)
  92     {
  93         const word& fieldName = fieldNames[i];
  94
  95         Info<< "Subsetting field " << fieldName << endl;
  96
  97         GeometricField<Type, fvsPatchField, surfaceMesh> fld
  98         (
  99             IOobject
 100             (
 101                 fieldName,
 102                 baseMesh.time().timeName(),
 103                 baseMesh,
 104                 IOobject::MUST_READ,
 105                 IOobject::NO_WRITE
 106             ),
 107             baseMesh
 108         );
 109
 110         subFields.set(i, subsetter.interpolate(fld));
 111     }
 112 }
 113
 114
 115 template<class Type>
 116 void subsetPointFields
 117 (
 118     const fvMeshSubset& subsetter,
 119     const pointMesh& pMesh,
 120     const wordList& fieldNames,
 121     PtrList<GeometricField<Type, pointPatchField, pointMesh> >& subFields
 122 )
 123 {
 124     const fvMesh& baseMesh = subsetter.baseMesh();
 125
 126     forAll(fieldNames, i)
 127     {
 128         const word& fieldName = fieldNames[i];
 129
 130         Info<< "Subsetting field " << fieldName << endl;
 131
 132         GeometricField<Type, pointPatchField, pointMesh> fld
 133         (
 134             IOobject
 135             (
 136                 fieldName,
 137                 baseMesh.time().timeName(),
 138                 baseMesh,
 139                 IOobject::MUST_READ,
 140                 IOobject::NO_WRITE
 141             ),
 142             pMesh
 143         );
 144
 145         subFields.set(i, subsetter.interpolate(fld));
 146     }
 147 }
 148
 149
 150 // Main program:
 151
 152 int main(int argc, char *argv[])
 153 {
 154     argList::validArgs.append("set");
 155     argList::validOptions.insert("patch", "patch name");
 156     argList::validOptions.insert("overwrite", "");
 157
 158 #   include "setRootCase.H"
 159 #   include "createTime.H"
 160 #   include "createMesh.H"
 161
 162     word setName(args.additionalArgs()[0]);
 163     bool overwrite = args.options().found("overwrite");
 164
 165
 166     Info<< "Reading cell set from " << setName << endl << endl;
 167
 168     // Create mesh subsetting engine
 169     fvMeshSubset subsetter(mesh);
 170
 171     label patchI = -1;
 172
 173     if (args.options().found("patch"))
 174     {
 175         word patchName(args.options()["patch"]);
 176
 177         patchI = mesh.boundaryMesh().findPatchID(patchName);
 178
 179         if (patchI == -1)
 180         {
 181             FatalErrorIn(args.executable()) << "Illegal patch " << patchName
 182                 << nl << "Valid patches are " << mesh.boundaryMesh().names()
 183                 << exit(FatalError);
 184         }
 185
 186         Info<< "Adding exposed internal faces to patch " << patchName << endl
 187             << endl;
 188     }
 189     else
 190     {
 191         Info<< "Adding exposed internal faces to a patch called"
 192             << " \"oldInternalFaces\" (created if nessecary)" << endl
 193             << endl;
 194     }
 195
 196
 197     cellSet currentSet(mesh, setName);
 198
 199     subsetter.setLargeCellSubset(currentSet, patchI, true);
 200
 201     IOobjectList objects(mesh, runTime.timeName());
 202
 203
 204     // Read vol fields and subset
 205     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 206
 207     wordList scalarNames(objects.names(volScalarField::typeName));
 208     PtrList<volScalarField> scalarFlds(scalarNames.size());
 209     subsetVolFields(subsetter, scalarNames, scalarFlds);
 210
 211     wordList vectorNames(objects.names(volVectorField::typeName));
 212     PtrList<volVectorField> vectorFlds(vectorNames.size());
 213     subsetVolFields(subsetter, vectorNames, vectorFlds);
 214
 215     wordList sphericalTensorNames
 216     (
 217         objects.names(volSphericalTensorField::typeName)
 218     );
 219     PtrList<volSphericalTensorField> sphericalTensorFlds
 220     (
 221         sphericalTensorNames.size()
 222     );
 223     subsetVolFields(subsetter, sphericalTensorNames, sphericalTensorFlds);
 224
 225     wordList symmTensorNames(objects.names(volSymmTensorField::typeName));
 226     PtrList<volSymmTensorField> symmTensorFlds(symmTensorNames.size());
 227     subsetVolFields(subsetter, symmTensorNames, symmTensorFlds);
 228
 229     wordList tensorNames(objects.names(volTensorField::typeName));
 230     PtrList<volTensorField> tensorFlds(tensorNames.size());
 231     subsetVolFields(subsetter, tensorNames, tensorFlds);
 232
 233
 234     // Read surface fields and subset
 235     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 236
 237     wordList surfScalarNames(objects.names(surfaceScalarField::typeName));
 238     PtrList<surfaceScalarField> surfScalarFlds(surfScalarNames.size());
 239     subsetSurfaceFields(subsetter, surfScalarNames, surfScalarFlds);
 240
 241     wordList surfVectorNames(objects.names(surfaceVectorField::typeName));
 242     PtrList<surfaceVectorField> surfVectorFlds(surfVectorNames.size());
 243     subsetSurfaceFields(subsetter, surfVectorNames, surfVectorFlds);
 244
 245     wordList surfSphericalTensorNames
 246     (
 247         objects.names(surfaceSphericalTensorField::typeName)
 248     );
 249     PtrList<surfaceSphericalTensorField> surfSphericalTensorFlds
 250     (
 251         surfSphericalTensorNames.size()
 252     );
 253     subsetSurfaceFields
 254     (
 255         subsetter,
 256         surfSphericalTensorNames,
 257         surfSphericalTensorFlds
 258     );
 259
 260     wordList surfSymmTensorNames
 261     (
 262         objects.names(surfaceSymmTensorField::typeName)
 263     );
 264     PtrList<surfaceSymmTensorField> surfSymmTensorFlds
 265     (
 266         surfSymmTensorNames.size()
 267     );
 268     subsetSurfaceFields(subsetter, surfSymmTensorNames, surfSymmTensorFlds);
 269
 270     wordList surfTensorNames(objects.names(surfaceTensorField::typeName));
 271     PtrList<surfaceTensorField> surfTensorFlds(surfTensorNames.size());
 272     subsetSurfaceFields(subsetter, surfTensorNames, surfTensorFlds);
 273
 274
 275     // Read point fields and subset
 276     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 277
 278     pointMesh pMesh(mesh);
 279
 280     wordList pointScalarNames(objects.names(pointScalarField::typeName));
 281     PtrList<pointScalarField> pointScalarFlds(pointScalarNames.size());
 282     subsetPointFields(subsetter, pMesh, pointScalarNames, pointScalarFlds);
 283
 284     wordList pointVectorNames(objects.names(pointVectorField::typeName));
 285     PtrList<pointVectorField> pointVectorFlds(pointVectorNames.size());
 286     subsetPointFields(subsetter, pMesh, pointVectorNames, pointVectorFlds);
 287
 288     wordList pointSphericalTensorNames
 289     (
 290         objects.names(pointSphericalTensorField::typeName)
 291     );
 292     PtrList<pointSphericalTensorField> pointSphericalTensorFlds
 293     (
 294         pointSphericalTensorNames.size()
 295     );
 296     subsetPointFields
 297     (
 298         subsetter,
 299         pMesh,
 300         pointSphericalTensorNames,
 301         pointSphericalTensorFlds
 302     );
 303
 304     wordList pointSymmTensorNames
 305     (
 306         objects.names(pointSymmTensorField::typeName)
 307     );
 308     PtrList<pointSymmTensorField> pointSymmTensorFlds
 309     (
 310         pointSymmTensorNames.size()
 311     );
 312     subsetPointFields
 313     (
 314         subsetter,
 315         pMesh,
 316         pointSymmTensorNames,
 317         pointSymmTensorFlds
 318     );
 319
 320     wordList pointTensorNames(objects.names(pointTensorField::typeName));
 321     PtrList<pointTensorField> pointTensorFlds(pointTensorNames.size());
 322     subsetPointFields(subsetter, pMesh, pointTensorNames, pointTensorFlds);
 323
 324
 325
 326     // Write mesh and fields to new time
 327     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 328
 329     if (!overwrite)
 330     {
 331         runTime++;
 332     }
 333
 334     Info<< "Writing subsetted mesh and fields to time " << runTime.value()
 335         << endl;
 336     subsetter.subMesh().write();
 337
 338
 339     // Subsetting adds 'subset' prefix. Rename field to be like original.
 340     forAll(scalarFlds, i)
 341     {
 342         scalarFlds[i].rename(scalarNames[i]);
 343
 344         scalarFlds[i].write();
 345     }
 346     forAll(vectorFlds, i)
 347     {
 348         vectorFlds[i].rename(vectorNames[i]);
 349
 350         vectorFlds[i].write();
 351     }
 352     forAll(sphericalTensorFlds, i)
 353     {
 354         sphericalTensorFlds[i].rename(sphericalTensorNames[i]);
 355
 356         sphericalTensorFlds[i].write();
 357     }
 358     forAll(symmTensorFlds, i)
 359     {
 360         symmTensorFlds[i].rename(symmTensorNames[i]);
 361
 362         symmTensorFlds[i].write();
 363     }
 364     forAll(tensorFlds, i)
 365     {
 366         tensorFlds[i].rename(tensorNames[i]);
 367
 368         tensorFlds[i].write();
 369     }
 370
 371     // Surface ones.
 372     forAll(surfScalarFlds, i)
 373     {
 374         surfScalarFlds[i].rename(surfScalarNames[i]);
 375
 376         surfScalarFlds[i].write();
 377     }
 378     forAll(surfVectorFlds, i)
 379     {
 380         surfVectorFlds[i].rename(surfVectorNames[i]);
 381
 382         surfVectorFlds[i].write();
 383     }
 384     forAll(surfSphericalTensorFlds, i)
 385     {
 386         surfSphericalTensorFlds[i].rename(surfSphericalTensorNames[i]);
 387
 388         surfSphericalTensorFlds[i].write();
 389     }
 390     forAll(surfSymmTensorFlds, i)
 391     {
 392         surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
 393
 394         surfSymmTensorFlds[i].write();
 395     }
 396     forAll(surfTensorNames, i)
 397     {
 398         surfTensorFlds[i].rename(surfTensorNames[i]);
 399
 400         surfTensorFlds[i].write();
 401     }
 402
 403     // Point ones
 404     forAll(pointScalarFlds, i)
 405     {
 406         pointScalarFlds[i].rename(pointScalarNames[i]);
 407
 408         pointScalarFlds[i].write();
 409     }
 410     forAll(pointVectorFlds, i)
 411     {
 412         pointVectorFlds[i].rename(pointVectorNames[i]);
 413
 414         pointVectorFlds[i].write();
 415     }
 416     forAll(pointSphericalTensorFlds, i)
 417     {
 418         pointSphericalTensorFlds[i].rename(pointSphericalTensorNames[i]);
 419
 420         pointSphericalTensorFlds[i].write();
 421     }
 422     forAll(pointSymmTensorFlds, i)
 423     {
 424         pointSymmTensorFlds[i].rename(pointSymmTensorNames[i]);
 425
 426         pointSymmTensorFlds[i].write();
 427     }
 428     forAll(pointTensorNames, i)
 429     {
 430         pointTensorFlds[i].rename(pointTensorNames[i]);
 431
 432         pointTensorFlds[i].write();
 433     }
 434
 435
 436     Info << nl << "End" << endl;
 437
 438     return 0;
 439 }
 440
 441
 442 // ************************************************************************* //