applications/utilities/parallelProcessing/reconstructParMesh/fvFieldReconstructorReconstructFields.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "fvFieldReconstructor.H"
  27 #include "foamTime.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                                 procField.boundaryField()[patchI].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 =
 161                             mesh_.boundaryMesh().whichPatch(curF);
 162
 163                         if (!patchFields(curBPatch))
 164                         {
 165                             patchFields.set
 166                             (
 167                                 curBPatch,
 168                                 fvPatchField<Type>::New
 169                                 (
 170                                     mesh_.boundary()[curBPatch].type(),
 171                                     mesh_.boundary()[curBPatch],
 172                                     DimensionedField<Type, volMesh>::null()
 173                                 )
 174                             );
 175                         }
 176
 177                         // add the face
 178                         label curPatchFace =
 179                             mesh_.boundaryMesh()
 180                                 [curBPatch].whichFace(curF);
 181
 182                         patchFields[curBPatch][curPatchFace] =
 183                             curProcPatch[faceI];
 184                     }
 185                 }
 186             }
 187         }
 188     }
 189
 190     forAll (mesh_.boundary(), patchI)
 191     {
 192         // add empty patches
 193         if
 194         (
 195             isType<emptyFvPatch>(mesh_.boundary()[patchI])
 196          && !patchFields(patchI)
 197         )
 198         {
 199             patchFields.set
 200             (
 201                 patchI,
 202                 fvPatchField<Type>::New
 203                 (
 204                     emptyFvPatchField<Type>::typeName,
 205                     mesh_.boundary()[patchI],
 206                     DimensionedField<Type, volMesh>::null()
 207                 )
 208             );
 209         }
 210     }
 211
 212
 213     // Now construct and write the field
 214     // setting the internalField and patchFields
 215     return tmp<GeometricField<Type, fvPatchField, volMesh> >
 216     (
 217         new GeometricField<Type, fvPatchField, volMesh>
 218         (
 219             IOobject
 220             (
 221                 fieldIoObject.name(),
 222                 mesh_.time().timeName(),
 223                 mesh_,
 224                 IOobject::NO_READ,
 225                 IOobject::NO_WRITE
 226             ),
 227             mesh_,
 228             procFields[0].dimensions(),
 229             internalField,
 230             patchFields
 231         )
 232     );
 233 }
 234
 235
 236 template<class Type>
 237 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 238 Foam::fvFieldReconstructor::reconstructFvSurfaceField
 239 (
 240     const IOobject& fieldIoObject
 241 )
 242 {
 243     // Read the field for all the processors
 244     PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
 245     (
 246         procMeshes_.size()
 247     );
 248
 249     forAll (procMeshes_, procI)
 250     {
 251         procFields.set
 252         (
 253             procI,
 254             new GeometricField<Type, fvsPatchField, surfaceMesh>
 255             (
 256                 IOobject
 257                 (
 258                     fieldIoObject.name(),
 259                     procMeshes_[procI].time().timeName(),
 260                     procMeshes_[procI],
 261                     IOobject::MUST_READ,
 262                     IOobject::NO_WRITE
 263                 ),
 264                 procMeshes_[procI]
 265             )
 266         );
 267     }
 268
 269
 270     // Create the internalField
 271     Field<Type> internalField(mesh_.nInternalFaces());
 272
 273     // Create the patch fields
 274     PtrList<fvsPatchField<Type> > patchFields(mesh_.boundary().size());
 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                 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                                 procField.boundaryField()[patchI].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             isType<emptyFvPatch>(mesh_.boundary()[patchI])
 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
 490 // Reconstruct and write all/selected surface fields
 491 template<class Type>
 492 void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
 493 (
 494     const IOobjectList& objects,
 495     const HashSet<word>& selectedFields
 496 )
 497 {
 498     const word& fieldClassName =
 499         GeometricField<Type, fvsPatchField, surfaceMesh>::typeName;
 500
 501     IOobjectList fields = objects.lookupClass(fieldClassName);
 502
 503     if (fields.size())
 504     {
 505         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 506
 507         forAllConstIter(IOobjectList, fields, fieldIter)
 508         {
 509             if
 510             (
 511                 !selectedFields.size()
 512              || selectedFields.found(fieldIter()->name())
 513             )
 514             {
 515                 Info<< "        " << fieldIter()->name() << endl;
 516
 517                 reconstructFvSurfaceField<Type>(*fieldIter())().write();
 518             }
 519         }
 520
 521         Info<< endl;
 522     }
 523 }
 524
 525
 526 // ************************************************************************* //