src/dynamicMesh/dynamicFvMesh/fvMeshAdder/fvMeshAdderTemplates.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
  28 #include "volFields.H"
  29 #include "surfaceFields.H"
  30 #include "emptyFvPatchField.H"
  31
  32 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  33
  34 template<class Type>
  35 void Foam::fvMeshAdder::map
  36 (
  37     const Field<Type>& oldFld,
  38     const labelList& oldToNew,
  39     Field<Type>& fld
  40 )
  41 {
  42     forAll(oldFld, cellI)
  43     {
  44         label newCellI = oldToNew[cellI];
  45
  46         if (newCellI >= 0 && newCellI < fld.size())
  47         {
  48             fld[newCellI] = oldFld[cellI];
  49         }
  50     }
  51 }
  52
  53
  54 template<class Type>
  55 void Foam::fvMeshAdder::MapVolField
  56 (
  57     const mapAddedPolyMesh& meshMap,
  58
  59     GeometricField<Type, fvPatchField, volMesh>& fld,
  60     const GeometricField<Type, fvPatchField, volMesh>& fldToAdd
  61 )
  62 {
  63     const fvMesh& mesh = fld.mesh();
  64
  65     // Internal field
  66     // ~~~~~~~~~~~~~~
  67
  68     {
  69         // Store old internal field
  70         Field<Type> oldInternalField(fld.internalField());
  71
  72         // Modify internal field
  73         Field<Type>& intFld = fld.internalField();
  74
  75         intFld.setSize(mesh.nCells());
  76
  77         map(oldInternalField, meshMap.oldCellMap(), intFld);
  78         map(fldToAdd.internalField(), meshMap.addedCellMap(), intFld);
  79     }
  80
  81
  82     // Patch fields from old mesh
  83     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
  84
  85     {
  86         const labelList& oldPatchMap = meshMap.oldPatchMap();
  87         const labelList& oldPatchStarts = meshMap.oldPatchStarts();
  88         const labelList& oldPatchSizes = meshMap.oldPatchSizes();
  89
  90         // Reorder old patches in order of new ones. Put removed patches at end.
  91
  92         label unusedPatchI = 0;
  93
  94         forAll(oldPatchMap, patchI)
  95         {
  96             label newPatchI = oldPatchMap[patchI];
  97
  98             if (newPatchI != -1)
  99             {
 100                 unusedPatchI++;
 101             }
 102         }
 103
 104         label nUsedPatches = unusedPatchI;
 105
 106         // Reorder list for patchFields
 107         labelList oldToNew(oldPatchMap.size());
 108
 109         forAll(oldPatchMap, patchI)
 110         {
 111             label newPatchI = oldPatchMap[patchI];
 112
 113             if (newPatchI != -1)
 114             {
 115                 oldToNew[patchI] = newPatchI;
 116             }
 117             else
 118             {
 119                 oldToNew[patchI] = unusedPatchI++;
 120             }
 121         }
 122
 123
 124         // Sort deleted ones last so is now in newPatch ordering
 125         fld.boundaryField().reorder(oldToNew);
 126         // Extend to covers all patches
 127         fld.boundaryField().setSize(mesh.boundaryMesh().size());
 128         // Delete unused patches
 129         for
 130         (
 131             label newPatchI = nUsedPatches;
 132             newPatchI < fld.boundaryField().size();
 133             newPatchI++
 134         )
 135         {
 136             fld.boundaryField().set(newPatchI, NULL);
 137         }
 138
 139
 140         // Map old values
 141         // ~~~~~~~~~~~~~~
 142
 143         forAll(oldPatchMap, patchI)
 144         {
 145             label newPatchI = oldPatchMap[patchI];
 146
 147             if (newPatchI != -1)
 148             {
 149                 labelList newToOld
 150                 (
 151                     calcPatchMap
 152                     (
 153                         oldPatchStarts[patchI],
 154                         oldPatchSizes[patchI],
 155                         meshMap.oldFaceMap(),
 156                         mesh.boundaryMesh()[newPatchI],
 157                         0                   // unmapped value
 158                     )
 159                 );
 160
 161                 directFvPatchFieldMapper patchMapper
 162                 (
 163                     oldPatchSizes[patchI],
 164                     newToOld
 165                 );
 166
 167
 168                 // Create new patchField with same type as existing one.
 169                 // Note:
 170                 // - boundaryField already in new order so access with newPatchI
 171                 // - fld.boundaryField()[newPatchI] both used for type and old
 172                 //   value
 173                 // - hope that field mapping allows aliasing since old and new
 174                 //   are same memory!
 175                 fld.boundaryField().set
 176                 (
 177                     newPatchI,
 178                     fvPatchField<Type>::New
 179                     (
 180                         fld.boundaryField()[newPatchI], // old field
 181                         mesh.boundary()[newPatchI],     // new fvPatch
 182                         fld.dimensionedInternalField(), // new internal field
 183                         patchMapper                     // mapper (new to old)
 184                     )
 185                 );
 186             }
 187         }
 188     }
 189
 190
 191
 192     // Patch fields from added mesh
 193     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 194
 195     {
 196         const labelList& addedPatchMap = meshMap.addedPatchMap();
 197
 198         // Add addedMesh patches
 199         forAll(addedPatchMap, patchI)
 200         {
 201             label newPatchI = addedPatchMap[patchI];
 202
 203             if (newPatchI != -1)
 204             {
 205                 const polyPatch& newPatch = mesh.boundaryMesh()[newPatchI];
 206                 const polyPatch& oldPatch =
 207                     fldToAdd.mesh().boundaryMesh()[patchI];
 208
 209                 if (!fld.boundaryField()(newPatchI))
 210                 {
 211                     // First occurrence of newPatchI. Map from existing
 212                     // patchField
 213
 214                     // From new patch faces to patch faces on added mesh.
 215                     labelList newToAdded
 216                     (
 217                         calcPatchMap
 218                         (
 219                             oldPatch.start(),
 220                             oldPatch.size(),
 221                             meshMap.addedFaceMap(),
 222                             newPatch,
 223                             0                       // unmapped values
 224                         )
 225                     );
 226
 227                     directFvPatchFieldMapper patchMapper
 228                     (
 229                         oldPatch.size(),
 230                         newToAdded
 231                     );
 232
 233                     fld.boundaryField().set
 234                     (
 235                         newPatchI,
 236                         fvPatchField<Type>::New
 237                         (
 238                             fldToAdd.boundaryField()[patchI], // added field
 239                             mesh.boundary()[newPatchI],       // new fvPatch
 240                             fld.dimensionedInternalField(),   // new int. field
 241                             patchMapper                       // mapper
 242                         )
 243                     );
 244                 }
 245                 else
 246                 {
 247                     // PatchField will have correct size already. Just slot in
 248                     // my elements.
 249
 250                     // From new patch faces to patch faces on added mesh. This
 251                     // time keep unmapped elements -1.
 252                     labelList newToAdded
 253                     (
 254                         calcPatchMap
 255                         (
 256                             oldPatch.start(),
 257                             oldPatch.size(),
 258                             meshMap.addedFaceMap(),
 259                             newPatch,
 260                             -1                      // unmapped values
 261                         )
 262                     );
 263
 264                     const fvPatchField<Type>& addedFld =
 265                         fldToAdd.boundaryField()[patchI];
 266
 267                     fvPatchField<Type>& newFld = fld.boundaryField()[newPatchI];
 268
 269                     forAll(newFld, i)
 270                     {
 271                         label oldFaceI = newToAdded[i];
 272
 273                         if (oldFaceI >= 0 && oldFaceI < addedFld.size())
 274                         {
 275                             newFld[i] = addedFld[oldFaceI];
 276                         }
 277                     }
 278                 }
 279             }
 280         }
 281     }
 282 }
 283
 284
 285 template<class Type>
 286 void Foam::fvMeshAdder::MapVolFields
 287 (
 288     const mapAddedPolyMesh& meshMap,
 289     const fvMesh& mesh,
 290     const fvMesh& meshToAdd
 291 )
 292 {
 293     HashTable<const GeometricField<Type, fvPatchField, volMesh>*> fields
 294     (
 295         mesh.objectRegistry::lookupClass
 296         <GeometricField<Type, fvPatchField, volMesh> >
 297         ()
 298     );
 299
 300     HashTable<const GeometricField<Type, fvPatchField, volMesh>*> fieldsToAdd
 301     (
 302         meshToAdd.objectRegistry::lookupClass
 303         <GeometricField<Type, fvPatchField, volMesh> >
 304         ()
 305     );
 306
 307     // It is necessary to enforce that all old-time fields are stored
 308     // before the mapping is performed.  Otherwise, if the
 309     // old-time-level field is mapped before the field itself, sizes
 310     // will not match.
 311
 312     for
 313     (
 314         typename HashTable<const GeometricField<Type, fvPatchField, volMesh>*>::
 315             iterator fieldIter = fields.begin();
 316         fieldIter != fields.end();
 317         ++fieldIter
 318     )
 319     {
 320         const_cast<GeometricField<Type, fvPatchField, volMesh>*>(fieldIter())
 321             ->storeOldTimes();
 322     }
 323
 324
 325     for
 326     (
 327         typename HashTable<const GeometricField<Type, fvPatchField, volMesh>*>::
 328             iterator fieldIter = fields.begin();
 329         fieldIter != fields.end();
 330         ++fieldIter
 331     )
 332     {
 333         GeometricField<Type, fvPatchField, volMesh>& fld =
 334             const_cast<GeometricField<Type, fvPatchField, volMesh>&>
 335             (
 336                 *fieldIter()
 337             );
 338
 339         if (fieldsToAdd.found(fld.name()))
 340         {
 341             Pout<< "Mapping field " << fld.name() << endl;
 342
 343             const GeometricField<Type, fvPatchField, volMesh>& fldToAdd =
 344                 *fieldsToAdd[fld.name()];
 345
 346             MapVolField<Type>(meshMap, fld, fldToAdd);
 347         }
 348         else
 349         {
 350             WarningIn("fvMeshAdder::MapVolFields")
 351                 << "Not mapping field " << fld.name()
 352                 << " since not present on mesh to add"
 353                 << endl;
 354         }
 355     }
 356 }
 357
 358
 359 template<class Type>
 360 void Foam::fvMeshAdder::MapSurfaceField
 361 (
 362     const mapAddedPolyMesh& meshMap,
 363
 364     GeometricField<Type, fvsPatchField, surfaceMesh>& fld,
 365     const GeometricField<Type, fvsPatchField, surfaceMesh>& fldToAdd
 366 )
 367 {
 368     const fvMesh& mesh = fld.mesh();
 369     const labelList& oldPatchStarts = meshMap.oldPatchStarts();
 370
 371     // Internal field
 372     // ~~~~~~~~~~~~~~
 373
 374     // Store old internal field
 375     {
 376         Field<Type> oldField(fld);
 377
 378         // Modify internal field
 379         Field<Type>& intFld = fld.internalField();
 380
 381         intFld.setSize(mesh.nInternalFaces());
 382
 383         map(oldField, meshMap.oldFaceMap(), intFld);
 384         map(fldToAdd, meshMap.addedFaceMap(), intFld);
 385
 386         // Faces that were boundary faces but are not anymore.
 387         // Use owner value (so lowest numbered cell, i.e. from 'old' not 'added'
 388         // mesh)
 389         forAll(fld.boundaryField(), patchI)
 390         {
 391             const fvsPatchField<Type>& pf = fld.boundaryField()[patchI];
 392
 393             label start = oldPatchStarts[patchI];
 394
 395             forAll(pf, i)
 396             {
 397                 label newFaceI = meshMap.oldFaceMap()[start + i];
 398
 399                 if (newFaceI >= 0 && newFaceI < mesh.nInternalFaces())
 400                 {
 401                     intFld[newFaceI] = pf[i];
 402                 }
 403             }
 404         }
 405     }
 406
 407
 408     // Patch fields from old mesh
 409     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 410
 411     {
 412         const labelList& oldPatchMap = meshMap.oldPatchMap();
 413         const labelList& oldPatchSizes = meshMap.oldPatchSizes();
 414
 415         // Reorder old patches in order of new ones. Put removed patches at end.
 416
 417         label unusedPatchI = 0;
 418
 419         forAll(oldPatchMap, patchI)
 420         {
 421             label newPatchI = oldPatchMap[patchI];
 422
 423             if (newPatchI != -1)
 424             {
 425                 unusedPatchI++;
 426             }
 427         }
 428
 429         label nUsedPatches = unusedPatchI;
 430
 431         // Reorder list for patchFields
 432         labelList oldToNew(oldPatchMap.size());
 433
 434         forAll(oldPatchMap, patchI)
 435         {
 436             label newPatchI = oldPatchMap[patchI];
 437
 438             if (newPatchI != -1)
 439             {
 440                 oldToNew[patchI] = newPatchI;
 441             }
 442             else
 443             {
 444                 oldToNew[patchI] = unusedPatchI++;
 445             }
 446         }
 447
 448
 449         // Sort deleted ones last so is now in newPatch ordering
 450         fld.boundaryField().reorder(oldToNew);
 451         // Extend to covers all patches
 452         fld.boundaryField().setSize(mesh.boundaryMesh().size());
 453         // Delete unused patches
 454         for
 455         (
 456             label newPatchI = nUsedPatches;
 457             newPatchI < fld.boundaryField().size();
 458             newPatchI++
 459         )
 460         {
 461             fld.boundaryField().set(newPatchI, NULL);
 462         }
 463
 464
 465         // Map old values
 466         // ~~~~~~~~~~~~~~
 467
 468         forAll(oldPatchMap, patchI)
 469         {
 470             label newPatchI = oldPatchMap[patchI];
 471
 472             if (newPatchI != -1)
 473             {
 474                 labelList newToOld
 475                 (
 476                     calcPatchMap
 477                     (
 478                         oldPatchStarts[patchI],
 479                         oldPatchSizes[patchI],
 480                         meshMap.oldFaceMap(),
 481                         mesh.boundaryMesh()[newPatchI],
 482                         0                   // unmapped value
 483                     )
 484                 );
 485
 486                 directFvPatchFieldMapper patchMapper
 487                 (
 488                     oldPatchSizes[patchI],
 489                     newToOld
 490                 );
 491
 492
 493                 // Create new patchField with same type as existing one.
 494                 // Note:
 495                 // - boundaryField already in new order so access with newPatchI
 496                 // - fld.boundaryField()[newPatchI] both used for type and old
 497                 //   value
 498                 // - hope that field mapping allows aliasing since old and new
 499                 //   are same memory!
 500                 fld.boundaryField().set
 501                 (
 502                     newPatchI,
 503                     fvsPatchField<Type>::New
 504                     (
 505                         fld.boundaryField()[newPatchI], // old field
 506                         mesh.boundary()[newPatchI],     // new fvPatch
 507                         fld.dimensionedInternalField(), // new internal field
 508                         patchMapper                     // mapper (new to old)
 509                     )
 510                 );
 511             }
 512         }
 513     }
 514
 515
 516
 517     // Patch fields from added mesh
 518     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 519
 520     {
 521         const labelList& addedPatchMap = meshMap.addedPatchMap();
 522
 523         // Add addedMesh patches
 524         forAll(addedPatchMap, patchI)
 525         {
 526             label newPatchI = addedPatchMap[patchI];
 527
 528             if (newPatchI != -1)
 529             {
 530                 const polyPatch& newPatch = mesh.boundaryMesh()[newPatchI];
 531                 const polyPatch& oldPatch =
 532                     fldToAdd.mesh().boundaryMesh()[patchI];
 533
 534                 if (!fld.boundaryField()(newPatchI))
 535                 {
 536                     // First occurrence of newPatchI. Map from existing
 537                     // patchField
 538
 539                     // From new patch faces to patch faces on added mesh.
 540                     labelList newToAdded
 541                     (
 542                         calcPatchMap
 543                         (
 544                             oldPatch.start(),
 545                             oldPatch.size(),
 546                             meshMap.addedFaceMap(),
 547                             newPatch,
 548                             0                       // unmapped values
 549                         )
 550                     );
 551
 552                     directFvPatchFieldMapper patchMapper
 553                     (
 554                         oldPatch.size(),
 555                         newToAdded
 556                     );
 557
 558                     fld.boundaryField().set
 559                     (
 560                         newPatchI,
 561                         fvsPatchField<Type>::New
 562                         (
 563                             fldToAdd.boundaryField()[patchI],// added field
 564                             mesh.boundary()[newPatchI],      // new fvPatch
 565                             fld.dimensionedInternalField(),  // new int. field
 566                             patchMapper                      // mapper
 567                         )
 568                     );
 569                 }
 570                 else
 571                 {
 572                     // PatchField will have correct size already. Just slot in
 573                     // my elements.
 574
 575                     // From new patch faces to patch faces on added mesh. This
 576                     // time keep unmapped elements -1.
 577                     labelList newToAdded
 578                     (
 579                         calcPatchMap
 580                         (
 581                             oldPatch.start(),
 582                             oldPatch.size(),
 583                             meshMap.addedFaceMap(),
 584                             newPatch,
 585                             -1                      // unmapped values
 586                         )
 587                     );
 588
 589                     const fvsPatchField<Type>& addedFld =
 590                         fldToAdd.boundaryField()[patchI];
 591
 592                     fvsPatchField<Type>& newFld =
 593                         fld.boundaryField()[newPatchI];
 594
 595                     forAll(newFld, i)
 596                     {
 597                         label oldFaceI = newToAdded[i];
 598
 599                         if (oldFaceI >= 0 && oldFaceI < addedFld.size())
 600                         {
 601                             newFld[i] = addedFld[oldFaceI];
 602                         }
 603                     }
 604                 }
 605             }
 606         }
 607     }
 608 }
 609
 610
 611 template<class Type>
 612 void Foam::fvMeshAdder::MapSurfaceFields
 613 (
 614     const mapAddedPolyMesh& meshMap,
 615     const fvMesh& mesh,
 616     const fvMesh& meshToAdd
 617 )
 618 {
 619     typedef GeometricField<Type, fvsPatchField, surfaceMesh> fldType;
 620
 621     HashTable<const fldType*> fields
 622     (
 623         mesh.objectRegistry::lookupClass<fldType>()
 624     );
 625
 626     HashTable<const fldType*> fieldsToAdd
 627     (
 628         meshToAdd.objectRegistry::lookupClass<fldType>()
 629     );
 630
 631     // It is necessary to enforce that all old-time fields are stored
 632     // before the mapping is performed.  Otherwise, if the
 633     // old-time-level field is mapped before the field itself, sizes
 634     // will not match.
 635
 636     for
 637     (
 638         typename HashTable<const fldType*>::
 639             iterator fieldIter = fields.begin();
 640         fieldIter != fields.end();
 641         ++fieldIter
 642     )
 643     {
 644         const_cast<fldType*>(fieldIter())
 645             ->storeOldTimes();
 646     }
 647
 648
 649     for
 650     (
 651         typename HashTable<const fldType*>::
 652             iterator fieldIter = fields.begin();
 653         fieldIter != fields.end();
 654         ++fieldIter
 655     )
 656     {
 657         fldType& fld = const_cast<fldType&>(*fieldIter());
 658
 659         if (fieldsToAdd.found(fld.name()))
 660         {
 661             Pout<< "Mapping field " << fld.name() << endl;
 662
 663             const fldType& fldToAdd = *fieldsToAdd[fld.name()];
 664
 665             MapSurfaceField<Type>(meshMap, fld, fldToAdd);
 666         }
 667         else
 668         {
 669             WarningIn("fvMeshAdder::MapSurfaceFields")
 670                 << "Not mapping field " << fld.name()
 671                 << " since not present on mesh to add"
 672                 << endl;
 673         }
 674     }
 675 }
 676
 677
 678 // ************************************************************************* //