applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C

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