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