applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-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 \*---------------------------------------------------------------------------*/
  25
  26 #include "fvFieldReconstructor.H"
  27 #include "Time.H"
  28 #include "PtrList.H"
  29 #include "fvPatchFields.H"
  30 #include "emptyFvPatch.H"
  31 #include "emptyFvPatchField.H"
  32 #include "emptyFvsPatchField.H"
  33
  34 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  35
  36 template<class Type>
  37 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
  38 Foam::fvFieldReconstructor::reconstructFvVolumeField
  39 (
  40     const IOobject& fieldIoObject
  41 )
  42 {
  43     // Read the field for all the processors
  44     PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
  45     (
  46         procMeshes_.size()
  47     );
  48
  49     forAll (procMeshes_, procI)
  50     {
  51         procFields.set
  52         (
  53             procI,
  54             new GeometricField<Type, fvPatchField, volMesh>
  55             (
  56                 IOobject
  57                 (
  58                     fieldIoObject.name(),
  59                     procMeshes_[procI].time().timeName(),
  60                     procMeshes_[procI],
  61                     IOobject::MUST_READ,
  62                     IOobject::NO_WRITE
  63                 ),
  64                 procMeshes_[procI]
  65             )
  66         );
  67     }
  68
  69
  70     // Create the internalField
  71     Field<Type> internalField(mesh_.nCells());
  72
  73     // Create the patch fields
  74     PtrList<fvPatchField<Type> > patchFields(mesh_.boundary().size());
  75
  76
  77     forAll (procMeshes_, procI)
  78     {
  79         const GeometricField<Type, fvPatchField, volMesh>& procField =
  80             procFields[procI];
  81
  82         // Set the cell values in the reconstructed field
  83         internalField.rmap
  84         (
  85             procField.internalField(),
  86             cellProcAddressing_[procI]
  87         );
  88
  89         // Set the boundary patch values in the reconstructed field
  90         forAll(boundaryProcAddressing_[procI], patchI)
  91         {
  92             // Get patch index of the original patch
  93             const label curBPatch = boundaryProcAddressing_[procI][patchI];
  94
  95             // Get addressing slice for this patch
  96             const labelList::subList cp =
  97                 procMeshes_[procI].boundary()[patchI].patchSlice
  98                 (
  99                     faceProcAddressing_[procI]
 100                 );
 101
 102             // check if the boundary patch is not a processor patch
 103             if (curBPatch >= 0)
 104             {
 105                 // Regular patch. Fast looping
 106
 107                 if (!patchFields(curBPatch))
 108                 {
 109                     patchFields.set
 110                     (
 111                         curBPatch,
 112                         fvPatchField<Type>::New
 113                         (
 114                             procField.boundaryField()[patchI],
 115                             mesh_.boundary()[curBPatch],
 116                             DimensionedField<Type, volMesh>::null(),
 117                             fvPatchFieldReconstructor
 118                             (
 119                                 mesh_.boundary()[curBPatch].size()
 120                             )
 121                         )
 122                     );
 123                 }
 124
 125                 const label curPatchStart =
 126                     mesh_.boundaryMesh()[curBPatch].start();
 127
 128                 labelList reverseAddressing(cp.size());
 129
 130                 forAll(cp, faceI)
 131                 {
 132                     // Subtract one to take into account offsets for
 133                     // face direction.
 134                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
 135                 }
 136
 137                 patchFields[curBPatch].rmap
 138                 (
 139                     procField.boundaryField()[patchI],
 140                     reverseAddressing
 141                 );
 142             }
 143             else
 144             {
 145                 const Field<Type>& curProcPatch =
 146                     procField.boundaryField()[patchI];
 147
 148                 // In processor patches, there's a mix of internal faces (some
 149                 // of them turned) and possible cyclics. Slow loop
 150                 forAll(cp, faceI)
 151                 {
 152                     // Subtract one to take into account offsets for
 153                     // face direction.
 154                     label curF = cp[faceI] - 1;
 155
 156                     // Is the face on the boundary?
 157                     if (curF >= mesh_.nInternalFaces())
 158                     {
 159                         label curBPatch = mesh_.boundaryMesh().whichPatch(curF);
 160
 161                         if (!patchFields(curBPatch))
 162                         {
 163                             patchFields.set
 164                             (
 165                                 curBPatch,
 166                                 fvPatchField<Type>::New
 167                                 (
 168                                     mesh_.boundary()[curBPatch].type(),
 169                                     mesh_.boundary()[curBPatch],
 170                                     DimensionedField<Type, volMesh>::null()
 171                                 )
 172                             );
 173                         }
 174
 175                         // add the face
 176                         label curPatchFace =
 177                             mesh_.boundaryMesh()
 178                                 [curBPatch].whichFace(curF);
 179
 180                         patchFields[curBPatch][curPatchFace] =
 181                             curProcPatch[faceI];
 182                     }
 183                 }
 184             }
 185         }
 186     }
 187
 188     forAll(mesh_.boundary(), patchI)
 189     {
 190         // add empty patches
 191         if
 192         (
 193             isType<emptyFvPatch>(mesh_.boundary()[patchI])
 194          && !patchFields(patchI)
 195         )
 196         {
 197             patchFields.set
 198             (
 199                 patchI,
 200                 fvPatchField<Type>::New
 201                 (
 202                     emptyFvPatchField<Type>::typeName,
 203                     mesh_.boundary()[patchI],
 204                     DimensionedField<Type, volMesh>::null()
 205                 )
 206             );
 207         }
 208     }
 209
 210
 211     // Now construct and write the field
 212     // setting the internalField and patchFields
 213     return tmp<GeometricField<Type, fvPatchField, volMesh> >
 214     (
 215         new GeometricField<Type, fvPatchField, volMesh>
 216         (
 217             IOobject
 218             (
 219                 fieldIoObject.name(),
 220                 mesh_.time().timeName(),
 221                 mesh_,
 222                 IOobject::NO_READ,
 223                 IOobject::NO_WRITE
 224             ),
 225             mesh_,
 226             procFields[0].dimensions(),
 227             internalField,
 228             patchFields
 229         )
 230     );
 231 }
 232
 233
 234 template<class Type>
 235 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 236 Foam::fvFieldReconstructor::reconstructFvSurfaceField
 237 (
 238     const IOobject& fieldIoObject
 239 )
 240 {
 241     // Read the field for all the processors
 242     PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
 243     (
 244         procMeshes_.size()
 245     );
 246
 247     forAll (procMeshes_, procI)
 248     {
 249         procFields.set
 250         (
 251             procI,
 252             new GeometricField<Type, fvsPatchField, surfaceMesh>
 253             (
 254                 IOobject
 255                 (
 256                     fieldIoObject.name(),
 257                     procMeshes_[procI].time().timeName(),
 258                     procMeshes_[procI],
 259                     IOobject::MUST_READ,
 260                     IOobject::NO_WRITE
 261                 ),
 262                 procMeshes_[procI]
 263             )
 264         );
 265     }
 266
 267
 268     // Create the internalField
 269     Field<Type> internalField(mesh_.nInternalFaces());
 270
 271     // Create the patch fields
 272     PtrList<fvsPatchField<Type> > patchFields(mesh_.boundary().size());
 273
 274
 275     forAll (procMeshes_, procI)
 276     {
 277         const GeometricField<Type, fvsPatchField, surfaceMesh>& procField =
 278             procFields[procI];
 279
 280         // Set the face values in the reconstructed field
 281
 282         // It is necessary to create a copy of the addressing array to
 283         // take care of the face direction offset trick.
 284         //
 285         {
 286             labelList curAddr(faceProcAddressing_[procI]);
 287
 288             forAll (curAddr, addrI)
 289             {
 290                 curAddr[addrI] -= 1;
 291             }
 292
 293             internalField.rmap
 294             (
 295                 procField.internalField(),
 296                 curAddr
 297             );
 298         }
 299
 300         // Set the boundary patch values in the reconstructed field
 301         forAll(boundaryProcAddressing_[procI], patchI)
 302         {
 303             // Get patch index of the original patch
 304             const label curBPatch = boundaryProcAddressing_[procI][patchI];
 305
 306             // Get addressing slice for this patch
 307             const labelList::subList cp =
 308                 procMeshes_[procI].boundary()[patchI].patchSlice
 309                 (
 310                     faceProcAddressing_[procI]
 311                 );
 312
 313             // check if the boundary patch is not a processor patch
 314             if (curBPatch >= 0)
 315             {
 316                 // Regular patch. Fast looping
 317
 318                 if (!patchFields(curBPatch))
 319                 {
 320                     patchFields.set
 321                     (
 322                         curBPatch,
 323                         fvsPatchField<Type>::New
 324                         (
 325                             procField.boundaryField()[patchI],
 326                             mesh_.boundary()[curBPatch],
 327                             DimensionedField<Type, surfaceMesh>::null(),
 328                             fvPatchFieldReconstructor
 329                             (
 330                                 mesh_.boundary()[curBPatch].size()
 331                             )
 332                         )
 333                     );
 334                 }
 335
 336                 const label curPatchStart =
 337                     mesh_.boundaryMesh()[curBPatch].start();
 338
 339                 labelList reverseAddressing(cp.size());
 340
 341                 forAll(cp, faceI)
 342                 {
 343                     // Subtract one to take into account offsets for
 344                     // face direction.
 345                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
 346                 }
 347
 348                 patchFields[curBPatch].rmap
 349                 (
 350                     procField.boundaryField()[patchI],
 351                     reverseAddressing
 352                 );
 353             }
 354             else
 355             {
 356                 const Field<Type>& curProcPatch =
 357                     procField.boundaryField()[patchI];
 358
 359                 // In processor patches, there's a mix of internal faces (some
 360                 // of them turned) and possible cyclics. Slow loop
 361                 forAll(cp, faceI)
 362                 {
 363                     label curF = cp[faceI] - 1;
 364
 365                     // Is the face turned the right side round
 366                     if (curF >= 0)
 367                     {
 368                         // Is the face on the boundary?
 369                         if (curF >= mesh_.nInternalFaces())
 370                         {
 371                             label curBPatch =
 372                                 mesh_.boundaryMesh().whichPatch(curF);
 373
 374                             if (!patchFields(curBPatch))
 375                             {
 376                                 patchFields.set
 377                                 (
 378                                     curBPatch,
 379                                     fvsPatchField<Type>::New
 380                                     (
 381                                         mesh_.boundary()[curBPatch].type(),
 382                                         mesh_.boundary()[curBPatch],
 383                                         DimensionedField<Type, surfaceMesh>
 384                                            ::null()
 385                                     )
 386                                 );
 387                             }
 388
 389                             // add the face
 390                             label curPatchFace =
 391                                 mesh_.boundaryMesh()
 392                                 [curBPatch].whichFace(curF);
 393
 394                             patchFields[curBPatch][curPatchFace] =
 395                                 curProcPatch[faceI];
 396                         }
 397                         else
 398                         {
 399                             // Internal face
 400                             internalField[curF] = curProcPatch[faceI];
 401                         }
 402                     }
 403                 }
 404             }
 405         }
 406     }
 407
 408     forAll(mesh_.boundary(), patchI)
 409     {
 410         // add empty patches
 411         if
 412         (
 413             isType<emptyFvPatch>(mesh_.boundary()[patchI])
 414          && !patchFields(patchI)
 415         )
 416         {
 417             patchFields.set
 418             (
 419                 patchI,
 420                 fvsPatchField<Type>::New
 421                 (
 422                     emptyFvsPatchField<Type>::typeName,
 423                     mesh_.boundary()[patchI],
 424                     DimensionedField<Type, surfaceMesh>::null()
 425                 )
 426             );
 427         }
 428     }
 429
 430
 431     // Now construct and write the field
 432     // setting the internalField and patchFields
 433     return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
 434     (
 435         new GeometricField<Type, fvsPatchField, surfaceMesh>
 436         (
 437             IOobject
 438             (
 439                 fieldIoObject.name(),
 440                 mesh_.time().timeName(),
 441                 mesh_,
 442                 IOobject::NO_READ,
 443                 IOobject::NO_WRITE
 444             ),
 445             mesh_,
 446             procFields[0].dimensions(),
 447             internalField,
 448             patchFields
 449         )
 450     );
 451 }
 452
 453
 454 // Reconstruct and write all/selected volume fields
 455 template<class Type>
 456 void Foam::fvFieldReconstructor::reconstructFvVolumeFields
 457 (
 458     const IOobjectList& objects,
 459     const HashSet<word>& selectedFields
 460 )
 461 {
 462     const word& fieldClassName =
 463         GeometricField<Type, fvPatchField, volMesh>::typeName;
 464
 465     IOobjectList fields = objects.lookupClass(fieldClassName);
 466
 467     if (fields.size())
 468     {
 469         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 470
 471         forAllConstIter(IOobjectList, fields, fieldIter)
 472         {
 473             if
 474             (
 475                 !selectedFields.size()
 476              || selectedFields.found(fieldIter()->name())
 477             )
 478             {
 479                 Info<< "        " << fieldIter()->name() << endl;
 480
 481                 reconstructFvVolumeField<Type>(*fieldIter())().write();
 482             }
 483         }
 484         Info<< endl;
 485     }
 486 }
 487
 488 // Reconstruct and write all/selected surface fields
 489 template<class Type>
 490 void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
 491 (
 492     const IOobjectList& objects,
 493     const HashSet<word>& selectedFields
 494 )
 495 {
 496     const word& fieldClassName =
 497         GeometricField<Type, fvsPatchField, surfaceMesh>::typeName;
 498
 499     IOobjectList fields = objects.lookupClass(fieldClassName);
 500
 501     if (fields.size())
 502     {
 503         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 504
 505         forAllConstIter(IOobjectList, fields, fieldIter)
 506         {
 507             if
 508             (
 509                 !selectedFields.size()
 510              || selectedFields.found(fieldIter()->name())
 511             )
 512             {
 513                 Info<< "        " << fieldIter()->name() << endl;
 514
 515                 reconstructFvSurfaceField<Type>(*fieldIter())().write();
 516             }
 517         }
 518         Info<< endl;
 519     }
 520 }
 521
 522
 523 // ************************************************************************* //