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