applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.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 "faFieldReconstructor.H"
  27 #include "foamTime.H"
  28 #include "PtrList.H"
  29 #include "faPatchFields.H"
  30 #include "emptyFaPatch.H"
  31 #include "emptyFaPatchField.H"
  32 #include "emptyFaePatchField.H"
  33
  34 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  35
  36 template<class Type>
  37 Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh> >
  38 Foam::faFieldReconstructor::reconstructFaAreaField
  39 (
  40     const IOobject& fieldIoObject
  41 )
  42 {
  43     // Read the field for all the processors
  44     PtrList<GeometricField<Type, faPatchField, areaMesh> > procFields
  45     (
  46         procMeshes_.size()
  47     );
  48
  49     forAll (procMeshes_, procI)
  50     {
  51         procFields.set
  52         (
  53             procI,
  54             new GeometricField<Type, faPatchField, areaMesh>
  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     // Create the internalField
  70     Field<Type> internalField(mesh_.nFaces());
  71
  72     // Create the patch fields
  73     PtrList<faPatchField<Type> > patchFields(mesh_.boundary().size());
  74
  75
  76     // Create global mesh patchs starts
  77
  78     labelList gStarts(mesh_.boundary().size(), -1);
  79
  80     if (mesh_.boundary().size() > 0)
  81     {
  82         gStarts[0] = mesh_.nInternalEdges();
  83     }
  84
  85     for(label i=1; i<mesh_.boundary().size(); i++)
  86     {
  87         gStarts[i] = gStarts[i-1] + mesh_.boundary()[i-1].labelList::size();
  88     }
  89
  90     forAll (procMeshes_, procI)
  91     {
  92         const GeometricField<Type, faPatchField, areaMesh>& procField =
  93             procFields[procI];
  94
  95         // Set the face values in the reconstructed field
  96         internalField.rmap
  97         (
  98             procField.internalField(),
  99             faceProcAddressing_[procI]
 100         );
 101
 102
 103
 104         // Set the boundary patch values in the reconstructed field
 105
 106         labelList starts(procMeshes_[procI].boundary().size(), -1);
 107
 108         if(procMeshes_[procI].boundary().size() > 0)
 109         {
 110             starts[0] = procMeshes_[procI].nInternalEdges();
 111         }
 112
 113         for(label i=1; i<procMeshes_[procI].boundary().size(); i++)
 114         {
 115             starts[i] =
 116                 starts[i-1]
 117               + procMeshes_[procI].boundary()[i-1].labelList::size();
 118         }
 119
 120         forAll(boundaryProcAddressing_[procI], patchI)
 121         {
 122             // Get patch index of the original patch
 123             const label curBPatch = boundaryProcAddressing_[procI][patchI];
 124
 125             // Get addressing slice for this patch
 126
 127 //             const labelList::subList cp =
 128 //                 procMeshes_[procI].boundary()[patchI].patchSlice
 129 //                 (
 130 //                     edgeProcAddressing_[procI]
 131 //                 );
 132
 133             const labelList::subList cp =
 134                 labelList::subList
 135                 (
 136                     edgeProcAddressing_[procI],
 137                     procMeshes_[procI].boundary()[patchI].size(),
 138                     starts[patchI]
 139                 );
 140
 141             // check if the boundary patch is not a processor patch
 142             if (curBPatch >= 0)
 143             {
 144                 // Regular patch. Fast looping
 145
 146                 if (!patchFields(curBPatch))
 147                 {
 148                     patchFields.set
 149                     (
 150                         curBPatch,
 151                         faPatchField<Type>::New
 152                         (
 153                             procField.boundaryField()[patchI],
 154                             mesh_.boundary()[curBPatch],
 155                             DimensionedField<Type, areaMesh>::null(),
 156                             faPatchFieldReconstructor
 157                             (
 158                                 mesh_.boundary()[curBPatch].size(),
 159                                 procField.boundaryField()[patchI].size()
 160                             )
 161                         )
 162                     );
 163                 }
 164
 165                 const label curPatchStart = gStarts[curBPatch];
 166 //                     mesh_.boundary()[curBPatch].start();
 167
 168                 labelList reverseAddressing(cp.size());
 169
 170                 forAll(cp, edgeI)
 171                 {
 172                     // Subtract one to take into account offsets for
 173                     // face direction.
 174 //                     reverseAddressing[edgeI] = cp[edgeI] - 1 - curPatchStart;
 175                     reverseAddressing[edgeI] = cp[edgeI] - curPatchStart;
 176                 }
 177
 178                 patchFields[curBPatch].rmap
 179                 (
 180                     procField.boundaryField()[patchI],
 181                     reverseAddressing
 182                 );
 183             }
 184             else
 185             {
 186                 const Field<Type>& curProcPatch =
 187                     procField.boundaryField()[patchI];
 188
 189                 // In processor patches, there's a mix of internal faces (some
 190                 // of them turned) and possible cyclics. Slow loop
 191                 forAll(cp, edgeI)
 192                 {
 193                     // Subtract one to take into account offsets for
 194                     // face direction.
 195 //                     label curE = cp[edgeI] - 1;
 196                     label curE = cp[edgeI];
 197
 198                     // Is the face on the boundary?
 199                     if (curE >= mesh_.nInternalEdges())
 200                     {
 201 //                         label curBPatch = mesh_.boundary().whichPatch(curE);
 202                         label curBPatch = -1;
 203
 204                         forAll (mesh_.boundary(), pI)
 205                         {
 206                             if
 207                             (
 208                                 curE >= gStarts[pI]
 209                              && curE <
 210                                 (
 211                                     gStarts[pI]
 212                                   + mesh_.boundary()[pI].labelList::size()
 213                                 )
 214                             )
 215                             {
 216                                 curBPatch = pI;
 217                             }
 218                         }
 219
 220                         if (!patchFields(curBPatch))
 221                         {
 222                             patchFields.set
 223                             (
 224                                 curBPatch,
 225                                 faPatchField<Type>::New
 226                                 (
 227                                     mesh_.boundary()[curBPatch].type(),
 228                                     mesh_.boundary()[curBPatch],
 229                                     DimensionedField<Type, areaMesh>::null()
 230                                 )
 231                             );
 232                         }
 233
 234                         // add the edge
 235 //                         label curPatchEdge =
 236 //                             mesh_.boundary()
 237 //                                 [curBPatch].whichEdge(curE);
 238
 239                         label curPatchEdge = curE - gStarts[curBPatch];
 240
 241                         patchFields[curBPatch][curPatchEdge] =
 242                             curProcPatch[edgeI];
 243                     }
 244                 }
 245             }
 246         }
 247     }
 248
 249     forAll(mesh_.boundary(), patchI)
 250     {
 251         // add empty patches
 252         if
 253         (
 254             typeid(mesh_.boundary()[patchI]) == typeid(emptyFaPatch)
 255          && !patchFields(patchI)
 256         )
 257         {
 258             patchFields.set
 259             (
 260                 patchI,
 261                 faPatchField<Type>::New
 262                 (
 263                     emptyFaPatchField<Type>::typeName,
 264                     mesh_.boundary()[patchI],
 265                     DimensionedField<Type, areaMesh>::null()
 266                 )
 267             );
 268         }
 269     }
 270
 271
 272     // Now construct and write the field
 273     // setting the internalField and patchFields
 274     return tmp<GeometricField<Type, faPatchField, areaMesh> >
 275     (
 276         new GeometricField<Type, faPatchField, areaMesh>
 277         (
 278             IOobject
 279             (
 280                 fieldIoObject.name(),
 281                 mesh_.time().timeName(),
 282                 mesh_(),
 283                 IOobject::NO_READ,
 284                 IOobject::NO_WRITE
 285             ),
 286             mesh_,
 287             procFields[0].dimensions(),
 288             internalField,
 289             patchFields
 290         )
 291     );
 292 }
 293
 294
 295 template<class Type>
 296 Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh> >
 297 Foam::faFieldReconstructor::reconstructFaEdgeField
 298 (
 299     const IOobject& fieldIoObject
 300 )
 301 {
 302     // Read the field for all the processors
 303     PtrList<GeometricField<Type, faePatchField, edgeMesh> > procFields
 304     (
 305         procMeshes_.size()
 306     );
 307
 308     forAll (procMeshes_, procI)
 309     {
 310         procFields.set
 311         (
 312             procI,
 313             new GeometricField<Type, faePatchField, edgeMesh>
 314             (
 315                 IOobject
 316                 (
 317                     fieldIoObject.name(),
 318                     procMeshes_[procI].time().timeName(),
 319                     procMeshes_[procI](),
 320                     IOobject::MUST_READ,
 321                     IOobject::NO_WRITE
 322                 ),
 323                 procMeshes_[procI]
 324             )
 325         );
 326     }
 327
 328
 329     // Create the internalField
 330     Field<Type> internalField(mesh_.nInternalEdges());
 331
 332     // Create the patch fields
 333     PtrList<faePatchField<Type> > patchFields(mesh_.boundary().size());
 334
 335
 336     labelList gStarts(mesh_.boundary().size(), -1);
 337
 338     if(mesh_.boundary().size() > 0)
 339     {
 340         gStarts[0] = mesh_.nInternalEdges();
 341     }
 342
 343     for(label i=1; i<mesh_.boundary().size(); i++)
 344     {
 345         gStarts[i] = gStarts[i-1] + mesh_.boundary()[i-1].labelList::size();
 346     }
 347
 348
 349     forAll (procMeshes_, procI)
 350     {
 351         const GeometricField<Type, faePatchField, edgeMesh>& procField =
 352             procFields[procI];
 353
 354         // Set the face values in the reconstructed field
 355
 356         // It is necessary to create a copy of the addressing array to
 357         // take care of the face direction offset trick.
 358         //
 359         {
 360             labelList curAddr(edgeProcAddressing_[procI]);
 361
 362 //             forAll (curAddr, addrI)
 363 //             {
 364 //                 curAddr[addrI] -= 1;
 365 //             }
 366
 367             internalField.rmap
 368             (
 369                 procField.internalField(),
 370                 curAddr
 371             );
 372         }
 373
 374         // Set the boundary patch values in the reconstructed field
 375
 376         labelList starts(procMeshes_[procI].boundary().size(), -1);
 377
 378         if(procMeshes_[procI].boundary().size() > 0)
 379         {
 380             starts[0] = procMeshes_[procI].nInternalEdges();
 381         }
 382
 383         for(label i=1; i<procMeshes_[procI].boundary().size(); i++)
 384         {
 385             starts[i] =
 386                 starts[i-1]
 387               + procMeshes_[procI].boundary()[i-1].labelList::size();
 388         }
 389
 390         forAll(boundaryProcAddressing_[procI], patchI)
 391         {
 392             // Get patch index of the original patch
 393             const label curBPatch = boundaryProcAddressing_[procI][patchI];
 394
 395             // Get addressing slice for this patch
 396
 397 //             const labelList::subList cp =
 398 //                 procMeshes_[procI].boundary()[patchI].patchSlice
 399 //                 (
 400 //                     faceProcAddressing_[procI]
 401 //                 );
 402
 403             const labelList::subList cp =
 404                 labelList::subList
 405                 (
 406                     edgeProcAddressing_[procI],
 407                     procMeshes_[procI].boundary()[patchI].size(),
 408                     starts[patchI]
 409                 );
 410
 411             // check if the boundary patch is not a processor patch
 412             if (curBPatch >= 0)
 413             {
 414                 // Regular patch. Fast looping
 415
 416                 if (!patchFields(curBPatch))
 417                 {
 418                     patchFields.set
 419                     (
 420                         curBPatch,
 421                         faePatchField<Type>::New
 422                         (
 423                             procField.boundaryField()[patchI],
 424                             mesh_.boundary()[curBPatch],
 425                             DimensionedField<Type, edgeMesh>::null(),
 426                             faPatchFieldReconstructor
 427                             (
 428                                 mesh_.boundary()[curBPatch].size(),
 429                                 procField.boundaryField()[patchI].size()
 430                             )
 431                         )
 432                     );
 433                 }
 434
 435                 const label curPatchStart = gStarts[curBPatch];
 436 //                     mesh_.boundary()[curBPatch].start();
 437
 438                 labelList reverseAddressing(cp.size());
 439
 440                 forAll(cp, edgeI)
 441                 {
 442                     // Subtract one to take into account offsets for
 443                     // face direction.
 444 //                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
 445                     reverseAddressing[edgeI] = cp[edgeI] - curPatchStart;
 446                 }
 447
 448                 patchFields[curBPatch].rmap
 449                 (
 450                     procField.boundaryField()[patchI],
 451                     reverseAddressing
 452                 );
 453             }
 454             else
 455             {
 456                 const Field<Type>& curProcPatch =
 457                     procField.boundaryField()[patchI];
 458
 459                 // In processor patches, there's a mix of internal faces (some
 460                 // of them turned) and possible cyclics. Slow loop
 461                 forAll(cp, edgeI)
 462                 {
 463 //                     label curF = cp[edgeI] - 1;
 464                     label curE = cp[edgeI];
 465
 466                     // Is the face turned the right side round
 467                     if (curE >= 0)
 468                     {
 469                         // Is the face on the boundary?
 470                         if (curE >= mesh_.nInternalEdges())
 471                         {
 472 //                             label curBPatch =
 473 //                                 mesh_.boundary().whichPatch(curF);
 474
 475                             label curBPatch = -1;
 476
 477                             forAll (mesh_.boundary(), pI)
 478                             {
 479                                 if
 480                                 (
 481                                     curE >= gStarts[pI]
 482                                  && curE <
 483                                     (
 484                                         gStarts[pI]
 485                                       + mesh_.boundary()[pI].labelList::size()
 486                                     )
 487                                 )
 488                                 {
 489                                     curBPatch = pI;
 490                                 }
 491                             }
 492
 493                             if (!patchFields(curBPatch))
 494                             {
 495                                 patchFields.set
 496                                 (
 497                                     curBPatch,
 498                                     faePatchField<Type>::New
 499                                     (
 500                                         mesh_.boundary()[curBPatch].type(),
 501                                         mesh_.boundary()[curBPatch],
 502                                         DimensionedField<Type, edgeMesh>
 503                                            ::null()
 504                                     )
 505                                 );
 506                             }
 507
 508                             // add the face
 509 //                             label curPatchFace =
 510 //                                 mesh_.boundary()
 511 //                                 [curBPatch].whichEdge(curF);
 512
 513                             label curPatchEdge = curE - gStarts[curBPatch];
 514
 515                             patchFields[curBPatch][curPatchEdge] =
 516                                 curProcPatch[edgeI];
 517                         }
 518                         else
 519                         {
 520                             // Internal face
 521                             internalField[curE] = curProcPatch[edgeI];
 522                         }
 523                     }
 524                 }
 525             }
 526         }
 527     }
 528
 529     forAll(mesh_.boundary(), patchI)
 530     {
 531         // add empty patches
 532         if
 533         (
 534             typeid(mesh_.boundary()[patchI]) == typeid(emptyFaPatch)
 535          && !patchFields(patchI)
 536         )
 537         {
 538             patchFields.set
 539             (
 540                 patchI,
 541                 faePatchField<Type>::New
 542                 (
 543                     emptyFaePatchField<Type>::typeName,
 544                     mesh_.boundary()[patchI],
 545                     DimensionedField<Type, edgeMesh>::null()
 546                 )
 547             );
 548         }
 549     }
 550
 551
 552     // Now construct and write the field
 553     // setting the internalField and patchFields
 554     return tmp<GeometricField<Type, faePatchField, edgeMesh> >
 555     (
 556         new GeometricField<Type, faePatchField, edgeMesh>
 557         (
 558             IOobject
 559             (
 560                 fieldIoObject.name(),
 561                 mesh_.time().timeName(),
 562                 mesh_(),
 563                 IOobject::NO_READ,
 564                 IOobject::NO_WRITE
 565             ),
 566             mesh_,
 567             procFields[0].dimensions(),
 568             internalField,
 569             patchFields
 570         )
 571     );
 572 }
 573
 574
 575 // Reconstruct and write all area fields
 576 template<class Type>
 577 void Foam::faFieldReconstructor::reconstructFaAreaFields
 578 (
 579     const IOobjectList& objects
 580 )
 581 {
 582     const word& fieldClassName =
 583         GeometricField<Type, faPatchField, areaMesh>::typeName;
 584
 585     IOobjectList fields = objects.lookupClass(fieldClassName);
 586
 587     if (fields.size())
 588     {
 589         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 590
 591         for
 592         (
 593             IOobjectList::const_iterator fieldIter = fields.begin();
 594             fieldIter != fields.end();
 595             ++fieldIter
 596         )
 597         {
 598             Info << "        " << fieldIter()->name() << endl;
 599
 600             reconstructFaAreaField<Type>(*fieldIter())().write();
 601         }
 602
 603         Info<< endl;
 604     }
 605 }
 606
 607 // Reconstruct and write all edge fields
 608 template<class Type>
 609 void Foam::faFieldReconstructor::reconstructFaEdgeFields
 610 (
 611     const IOobjectList& objects
 612 )
 613 {
 614     const word& fieldClassName =
 615         GeometricField<Type, faePatchField, edgeMesh>::typeName;
 616
 617     IOobjectList fields = objects.lookupClass(fieldClassName);
 618
 619     if (fields.size())
 620     {
 621         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 622
 623         for
 624         (
 625             IOobjectList::const_iterator fieldIter = fields.begin();
 626             fieldIter != fields.end();
 627             ++fieldIter
 628         )
 629         {
 630             Info<< "        " << fieldIter()->name() << endl;
 631
 632             reconstructFaEdgeField<Type>(*fieldIter())().write();
 633         }
 634
 635         Info<< endl;
 636     }
 637 }
 638
 639
 640 // ************************************************************************* //