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