applications/utilities/preProcessing/setFields/setFields.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     Set values on a selected set of cells/patchfaces through a dictionary.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "argList.H"
  30 #include "timeSelector.H"
  31 #include "Time.H"
  32 #include "fvMesh.H"
  33 #include "topoSetSource.H"
  34 #include "cellSet.H"
  35 #include "faceSet.H"
  36 #include "volFields.H"
  37
  38 using namespace Foam;
  39
  40 template<class Type>
  41 bool setCellFieldType
  42 (
  43     const word& fieldTypeDesc,
  44     const fvMesh& mesh,
  45     const labelList& selectedCells,
  46     Istream& fieldValueStream
  47 )
  48 {
  49     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
  50
  51     if (fieldTypeDesc != fieldType::typeName + "Value")
  52     {
  53         return false;
  54     }
  55
  56     word fieldName(fieldValueStream);
  57
  58     IOobject fieldHeader
  59     (
  60         fieldName,
  61         mesh.time().timeName(),
  62         mesh,
  63         IOobject::MUST_READ
  64     );
  65
  66     // Check field exists
  67     if (fieldHeader.headerOk())
  68     {
  69         Info<< "    Setting internal values of "
  70             << fieldHeader.headerClassName()
  71             << " " << fieldName << endl;
  72
  73         fieldType field(fieldHeader, mesh);
  74
  75         const Type& value = pTraits<Type>(fieldValueStream);
  76
  77         if (selectedCells.size() == field.size())
  78         {
  79             field.internalField() = value;
  80         }
  81         else
  82         {
  83             forAll(selectedCells, celli)
  84             {
  85                 field[selectedCells[celli]] = value;
  86             }
  87         }
  88
  89         forAll(field.boundaryField(), patchi)
  90         {
  91             field.boundaryField()[patchi] =
  92                 field.boundaryField()[patchi].patchInternalField();
  93         }
  94
  95         field.write();
  96     }
  97     else
  98     {
  99         WarningIn
 100         (
 101             "void setCellFieldType"
 102             "(const fvMesh& mesh, const labelList& selectedCells,"
 103             "Istream& fieldValueStream)"
 104         ) << "Field " << fieldName << " not found" << endl;
 105     }
 106
 107     return true;
 108 }
 109
 110
 111 class setCellField
 112 {
 113
 114 public:
 115
 116     setCellField()
 117     {}
 118
 119     autoPtr<setCellField> clone() const
 120     {
 121         return autoPtr<setCellField>(new setCellField());
 122     }
 123
 124     class iNew
 125     {
 126         const fvMesh& mesh_;
 127         const labelList& selectedCells_;
 128
 129     public:
 130
 131         iNew(const fvMesh& mesh, const labelList& selectedCells)
 132         :
 133             mesh_(mesh),
 134             selectedCells_(selectedCells)
 135         {}
 136
 137         autoPtr<setCellField> operator()(Istream& fieldValues) const
 138         {
 139             word fieldType(fieldValues);
 140
 141             if
 142             (
 143                !(
 144                     setCellFieldType<scalar>
 145                         (fieldType, mesh_, selectedCells_, fieldValues)
 146                  || setCellFieldType<vector>
 147                         (fieldType, mesh_, selectedCells_, fieldValues)
 148                  || setCellFieldType<sphericalTensor>
 149                         (fieldType, mesh_, selectedCells_, fieldValues)
 150                  || setCellFieldType<symmTensor>
 151                         (fieldType, mesh_, selectedCells_, fieldValues)
 152                  || setCellFieldType<tensor>
 153                         (fieldType, mesh_, selectedCells_, fieldValues)
 154                 )
 155             )
 156             {
 157                 WarningIn("setCellField::iNew::operator()(Istream& is)")
 158                     << "field type " << fieldType << " not currently supported"
 159                     << endl;
 160             }
 161
 162             return autoPtr<setCellField>(new setCellField());
 163         }
 164     };
 165 };
 166
 167
 168 template<class Type>
 169 bool setFaceFieldType
 170 (
 171     const word& fieldTypeDesc,
 172     const fvMesh& mesh,
 173     const labelList& selectedFaces,
 174     Istream& fieldValueStream
 175 )
 176 {
 177     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 178
 179     if (fieldTypeDesc != fieldType::typeName + "Value")
 180     {
 181         return false;
 182     }
 183
 184     word fieldName(fieldValueStream);
 185
 186     IOobject fieldHeader
 187     (
 188         fieldName,
 189         mesh.time().timeName(),
 190         mesh,
 191         IOobject::MUST_READ
 192     );
 193
 194     // Check field exists
 195     if (fieldHeader.headerOk())
 196     {
 197         Info<< "    Setting patchField values of "
 198             << fieldHeader.headerClassName()
 199             << " " << fieldName << endl;
 200
 201         fieldType field(fieldHeader, mesh);
 202
 203         const Type& value = pTraits<Type>(fieldValueStream);
 204
 205         // Create flat list of selected faces and their value.
 206         Field<Type> allBoundaryValues(mesh.nFaces()-mesh.nInternalFaces());
 207         forAll(field.boundaryField(), patchi)
 208         {
 209             SubField<Type>
 210             (
 211                 allBoundaryValues,
 212                 field.boundaryField()[patchi].size(),
 213                 field.boundaryField()[patchi].patch().start()
 214               - mesh.nInternalFaces()
 215             ).assign(field.boundaryField()[patchi]);
 216         }
 217
 218         // Override
 219         labelList nChanged(field.boundaryField().size(), 0);
 220         forAll(selectedFaces, i)
 221         {
 222             label facei = selectedFaces[i];
 223             if (mesh.isInternalFace(facei))
 224             {
 225                 WarningIn("setFaceFieldType(..)")
 226                     << "Ignoring internal face " << facei << endl;
 227             }
 228             else
 229             {
 230                 label bFaceI = facei-mesh.nInternalFaces();
 231                 allBoundaryValues[bFaceI] = value;
 232                 nChanged[mesh.boundaryMesh().patchID()[bFaceI]]++;
 233             }
 234         }
 235
 236         Pstream::listCombineGather(nChanged, plusEqOp<label>());
 237         Pstream::listCombineScatter(nChanged);
 238
 239         // Reassign.
 240         forAll(field.boundaryField(), patchi)
 241         {
 242             if (nChanged[patchi] > 0)
 243             {
 244                 Info<< "    On patch "
 245                     << field.boundaryField()[patchi].patch().name()
 246                     << " set " << nChanged[patchi] << " values" << endl;
 247                 field.boundaryField()[patchi] == SubField<Type>
 248                 (
 249                     allBoundaryValues,
 250                     field.boundaryField()[patchi].size(),
 251                     field.boundaryField()[patchi].patch().start()
 252                   - mesh.nInternalFaces()
 253                 );
 254             }
 255         }
 256
 257         field.write();
 258     }
 259     else
 260     {
 261         WarningIn
 262         (
 263             "void setFaceFieldType"
 264             "(const fvMesh& mesh, const labelList& selectedFaces,"
 265             "Istream& fieldValueStream)"
 266         ) << "Field " << fieldName << " not found" << endl;
 267     }
 268
 269     return true;
 270 }
 271
 272
 273 class setFaceField
 274 {
 275
 276 public:
 277
 278     setFaceField()
 279     {}
 280
 281     autoPtr<setFaceField> clone() const
 282     {
 283         return autoPtr<setFaceField>(new setFaceField());
 284     }
 285
 286     class iNew
 287     {
 288         const fvMesh& mesh_;
 289         const labelList& selectedFaces_;
 290
 291     public:
 292
 293         iNew(const fvMesh& mesh, const labelList& selectedFaces)
 294         :
 295             mesh_(mesh),
 296             selectedFaces_(selectedFaces)
 297         {}
 298
 299         autoPtr<setFaceField> operator()(Istream& fieldValues) const
 300         {
 301             word fieldType(fieldValues);
 302
 303             if
 304             (
 305                !(
 306                     setFaceFieldType<scalar>
 307                         (fieldType, mesh_, selectedFaces_, fieldValues)
 308                  || setFaceFieldType<vector>
 309                         (fieldType, mesh_, selectedFaces_, fieldValues)
 310                  || setFaceFieldType<sphericalTensor>
 311                         (fieldType, mesh_, selectedFaces_, fieldValues)
 312                  || setFaceFieldType<symmTensor>
 313                         (fieldType, mesh_, selectedFaces_, fieldValues)
 314                  || setFaceFieldType<tensor>
 315                         (fieldType, mesh_, selectedFaces_, fieldValues)
 316                 )
 317             )
 318             {
 319                 WarningIn("setFaceField::iNew::operator()(Istream& is)")
 320                     << "field type " << fieldType << " not currently supported"
 321                     << endl;
 322             }
 323
 324             return autoPtr<setFaceField>(new setFaceField());
 325         }
 326     };
 327 };
 328
 329
 330
 331 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 332
 333 int main(int argc, char *argv[])
 334 {
 335     timeSelector::addOptions();
 336
 337 #   include "setRootCase.H"
 338 #   include "createTime.H"
 339
 340     // Get times list
 341     instantList timeDirs = timeSelector::select0(runTime, args);
 342
 343 #   include "createMesh.H"
 344
 345     Info<< "Reading setFieldsDict\n" << endl;
 346
 347     IOdictionary setFieldsDict
 348     (
 349         IOobject
 350         (
 351             "setFieldsDict",
 352             runTime.system(),
 353             mesh,
 354             IOobject::MUST_READ_IF_MODIFIED,
 355             IOobject::NO_WRITE
 356         )
 357     );
 358
 359     if (setFieldsDict.found("defaultFieldValues"))
 360     {
 361         Info<< "Setting field default values" << endl;
 362         PtrList<setCellField> defaultFieldValues
 363         (
 364             setFieldsDict.lookup("defaultFieldValues"),
 365             setCellField::iNew(mesh, labelList(mesh.nCells()))
 366         );
 367         Info<< endl;
 368     }
 369
 370
 371     Info<< "Setting field region values" << endl;
 372
 373     PtrList<entry> regions(setFieldsDict.lookup("regions"));
 374
 375     forAll(regions, regionI)
 376     {
 377         const entry& region = regions[regionI];
 378
 379         autoPtr<topoSetSource> source =
 380             topoSetSource::New(region.keyword(), mesh, region.dict());
 381
 382         if (source().setType() == topoSetSource::CELLSETSOURCE)
 383         {
 384             cellSet selectedCellSet
 385             (
 386                 mesh,
 387                 "cellSet",
 388                 mesh.nCells()/10+1  // Reasonable size estimate.
 389             );
 390
 391             source->applyToSet
 392             (
 393                 topoSetSource::NEW,
 394                 selectedCellSet
 395             );
 396
 397             PtrList<setCellField> fieldValues
 398             (
 399                 region.dict().lookup("fieldValues"),
 400                 setCellField::iNew(mesh, selectedCellSet.toc())
 401             );
 402         }
 403         else if (source().setType() == topoSetSource::FACESETSOURCE)
 404         {
 405             faceSet selectedFaceSet
 406             (
 407                 mesh,
 408                 "faceSet",
 409                 (mesh.nFaces()-mesh.nInternalFaces())/10+1
 410             );
 411
 412             source->applyToSet
 413             (
 414                 topoSetSource::NEW,
 415                 selectedFaceSet
 416             );
 417
 418             PtrList<setFaceField> fieldValues
 419             (
 420                 region.dict().lookup("fieldValues"),
 421                 setFaceField::iNew(mesh, selectedFaceSet.toc())
 422             );
 423         }
 424     }
 425
 426
 427     Info<< "\nEnd\n" << endl;
 428
 429     return 0;
 430 }
 431
 432
 433 // ************************************************************************* //