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