src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   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(newToOld);
 162
 163
 164                 // Create new patchField with same type as existing one.
 165                 // Note:
 166                 // - boundaryField already in new order so access with newPatchI
 167                 // - fld.boundaryField()[newPatchI] both used for type and old
 168                 //   value
 169                 // - hope that field mapping allows aliasing since old and new
 170                 //   are same memory!
 171                 fld.boundaryField().set
 172                 (
 173                     newPatchI,
 174                     fvPatchField<Type>::New
 175                     (
 176                         fld.boundaryField()[newPatchI], // old field
 177                         mesh.boundary()[newPatchI],     // new fvPatch
 178                         fld.dimensionedInternalField(), // new internal field
 179                         patchMapper                     // mapper (new to old)
 180                     )
 181                 );
 182             }
 183         }
 184     }
 185
 186
 187
 188     // Patch fields from added mesh
 189     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 190
 191     {
 192         const labelList& addedPatchMap = meshMap.addedPatchMap();
 193
 194         // Add addedMesh patches
 195         forAll(addedPatchMap, patchI)
 196         {
 197             label newPatchI = addedPatchMap[patchI];
 198
 199             if (newPatchI != -1)
 200             {
 201                 const polyPatch& newPatch = mesh.boundaryMesh()[newPatchI];
 202                 const polyPatch& oldPatch =
 203                     fldToAdd.mesh().boundaryMesh()[patchI];
 204
 205                 if (!fld.boundaryField()(newPatchI))
 206                 {
 207                     // First occurrence of newPatchI. Map from existing
 208                     // patchField
 209
 210                     // From new patch faces to patch faces on added mesh.
 211                     labelList newToAdded
 212                     (
 213                         calcPatchMap
 214                         (
 215                             oldPatch.start(),
 216                             oldPatch.size(),
 217                             meshMap.addedFaceMap(),
 218                             newPatch,
 219                             0                       // unmapped values
 220                         )
 221                     );
 222
 223                     directFvPatchFieldMapper patchMapper(newToAdded);
 224
 225                     fld.boundaryField().set
 226                     (
 227                         newPatchI,
 228                         fvPatchField<Type>::New
 229                         (
 230                             fldToAdd.boundaryField()[patchI], // added field
 231                             mesh.boundary()[newPatchI],       // new fvPatch
 232                             fld.dimensionedInternalField(),   // new int. field
 233                             patchMapper                       // mapper
 234                         )
 235                     );
 236                 }
 237                 else
 238                 {
 239                     // PatchField will have correct size already. Just slot in
 240                     // my elements.
 241
 242                     // From new patch faces to patch faces on added mesh. This
 243                     // time keep unmapped elements -1.
 244                     labelList newToAdded
 245                     (
 246                         calcPatchMap
 247                         (
 248                             oldPatch.start(),
 249                             oldPatch.size(),
 250                             meshMap.addedFaceMap(),
 251                             newPatch,
 252                             -1                      // unmapped values
 253                         )
 254                     );
 255
 256                     const fvPatchField<Type>& addedFld =
 257                         fldToAdd.boundaryField()[patchI];
 258
 259                     fvPatchField<Type>& newFld = fld.boundaryField()[newPatchI];
 260
 261                     forAll(newFld, i)
 262                     {
 263                         label oldFaceI = newToAdded[i];
 264
 265                         if (oldFaceI >= 0 && oldFaceI < addedFld.size())
 266                         {
 267                             newFld[i] = addedFld[oldFaceI];
 268                         }
 269                     }
 270                 }
 271             }
 272         }
 273     }
 274 }
 275
 276
 277 template<class Type>
 278 void Foam::fvMeshAdder::MapVolFields
 279 (
 280     const mapAddedPolyMesh& meshMap,
 281     const fvMesh& mesh,
 282     const fvMesh& meshToAdd
 283 )
 284 {
 285     HashTable<const GeometricField<Type, fvPatchField, volMesh>*> fields
 286     (
 287         mesh.objectRegistry::lookupClass
 288         <GeometricField<Type, fvPatchField, volMesh> >
 289         ()
 290     );
 291
 292     HashTable<const GeometricField<Type, fvPatchField, volMesh>*> fieldsToAdd
 293     (
 294         meshToAdd.objectRegistry::lookupClass
 295         <GeometricField<Type, fvPatchField, volMesh> >
 296         ()
 297     );
 298
 299     // It is necessary to enforce that all old-time fields are stored
 300     // before the mapping is performed.  Otherwise, if the
 301     // old-time-level field is mapped before the field itself, sizes
 302     // will not match.
 303
 304     for
 305     (
 306         typename HashTable<const GeometricField<Type, fvPatchField, volMesh>*>::
 307             iterator fieldIter = fields.begin();
 308         fieldIter != fields.end();
 309         ++fieldIter
 310     )
 311     {
 312         const_cast<GeometricField<Type, fvPatchField, volMesh>*>(fieldIter())
 313             ->storeOldTimes();
 314     }
 315
 316
 317     for
 318     (
 319         typename HashTable<const GeometricField<Type, fvPatchField, volMesh>*>::
 320             iterator fieldIter = fields.begin();
 321         fieldIter != fields.end();
 322         ++fieldIter
 323     )
 324     {
 325         GeometricField<Type, fvPatchField, volMesh>& fld =
 326             const_cast<GeometricField<Type, fvPatchField, volMesh>&>
 327             (
 328                 *fieldIter()
 329             );
 330
 331         if (fieldsToAdd.found(fld.name()))
 332         {
 333             Pout<< "Mapping field " << fld.name() << endl;
 334
 335             const GeometricField<Type, fvPatchField, volMesh>& fldToAdd =
 336                 *fieldsToAdd[fld.name()];
 337
 338             MapVolField<Type>(meshMap, fld, fldToAdd);
 339         }
 340         else
 341         {
 342             WarningIn("fvMeshAdder::MapVolFields")
 343                 << "Not mapping field " << fld.name()
 344                 << " since not present on mesh to add"
 345                 << endl;
 346         }
 347     }
 348 }
 349
 350
 351 template<class Type>
 352 void Foam::fvMeshAdder::MapSurfaceField
 353 (
 354     const mapAddedPolyMesh& meshMap,
 355
 356     GeometricField<Type, fvsPatchField, surfaceMesh>& fld,
 357     const GeometricField<Type, fvsPatchField, surfaceMesh>& fldToAdd
 358 )
 359 {
 360     const fvMesh& mesh = fld.mesh();
 361     const labelList& oldPatchStarts = meshMap.oldPatchStarts();
 362
 363     // Internal field
 364     // ~~~~~~~~~~~~~~
 365
 366     // Store old internal field
 367     {
 368         Field<Type> oldField(fld);
 369
 370         // Modify internal field
 371         Field<Type>& intFld = fld.internalField();
 372
 373         intFld.setSize(mesh.nInternalFaces());
 374
 375         map(oldField, meshMap.oldFaceMap(), intFld);
 376         map(fldToAdd, meshMap.addedFaceMap(), intFld);
 377
 378         // Faces that were boundary faces but are not anymore.
 379         // Use owner value (so lowest numbered cell, i.e. from 'old' not 'added'
 380         // mesh)
 381         forAll(fld.boundaryField(), patchI)
 382         {
 383             const fvsPatchField<Type>& pf = fld.boundaryField()[patchI];
 384
 385             label start = oldPatchStarts[patchI];
 386
 387             forAll(pf, i)
 388             {
 389                 label newFaceI = meshMap.oldFaceMap()[start + i];
 390
 391                 if (newFaceI >= 0 && newFaceI < mesh.nInternalFaces())
 392                 {
 393                     intFld[newFaceI] = pf[i];
 394                 }
 395             }
 396         }
 397     }
 398
 399
 400     // Patch fields from old mesh
 401     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 402
 403     {
 404         const labelList& oldPatchMap = meshMap.oldPatchMap();
 405         const labelList& oldPatchSizes = meshMap.oldPatchSizes();
 406
 407         // Reorder old patches in order of new ones. Put removed patches at end.
 408
 409         label unusedPatchI = 0;
 410
 411         forAll(oldPatchMap, patchI)
 412         {
 413             label newPatchI = oldPatchMap[patchI];
 414
 415             if (newPatchI != -1)
 416             {
 417                 unusedPatchI++;
 418             }
 419         }
 420
 421         label nUsedPatches = unusedPatchI;
 422
 423         // Reorder list for patchFields
 424         labelList oldToNew(oldPatchMap.size());
 425
 426         forAll(oldPatchMap, patchI)
 427         {
 428             label newPatchI = oldPatchMap[patchI];
 429
 430             if (newPatchI != -1)
 431             {
 432                 oldToNew[patchI] = newPatchI;
 433             }
 434             else
 435             {
 436                 oldToNew[patchI] = unusedPatchI++;
 437             }
 438         }
 439
 440
 441         // Sort deleted ones last so is now in newPatch ordering
 442         fld.boundaryField().reorder(oldToNew);
 443         // Extend to covers all patches
 444         fld.boundaryField().setSize(mesh.boundaryMesh().size());
 445         // Delete unused patches
 446         for
 447         (
 448             label newPatchI = nUsedPatches;
 449             newPatchI < fld.boundaryField().size();
 450             newPatchI++
 451         )
 452         {
 453             fld.boundaryField().set(newPatchI, NULL);
 454         }
 455
 456
 457         // Map old values
 458         // ~~~~~~~~~~~~~~
 459
 460         forAll(oldPatchMap, patchI)
 461         {
 462             label newPatchI = oldPatchMap[patchI];
 463
 464             if (newPatchI != -1)
 465             {
 466                 labelList newToOld
 467                 (
 468                     calcPatchMap
 469                     (
 470                         oldPatchStarts[patchI],
 471                         oldPatchSizes[patchI],
 472                         meshMap.oldFaceMap(),
 473                         mesh.boundaryMesh()[newPatchI],
 474                         0                   // unmapped value
 475                     )
 476                 );
 477
 478                 directFvPatchFieldMapper patchMapper(newToOld);
 479
 480
 481                 // Create new patchField with same type as existing one.
 482                 // Note:
 483                 // - boundaryField already in new order so access with newPatchI
 484                 // - fld.boundaryField()[newPatchI] both used for type and old
 485                 //   value
 486                 // - hope that field mapping allows aliasing since old and new
 487                 //   are same memory!
 488                 fld.boundaryField().set
 489                 (
 490                     newPatchI,
 491                     fvsPatchField<Type>::New
 492                     (
 493                         fld.boundaryField()[newPatchI], // old field
 494                         mesh.boundary()[newPatchI],     // new fvPatch
 495                         fld.dimensionedInternalField(), // new internal field
 496                         patchMapper                     // mapper (new to old)
 497                     )
 498                 );
 499             }
 500         }
 501     }
 502
 503
 504
 505     // Patch fields from added mesh
 506     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 507
 508     {
 509         const labelList& addedPatchMap = meshMap.addedPatchMap();
 510
 511         // Add addedMesh patches
 512         forAll(addedPatchMap, patchI)
 513         {
 514             label newPatchI = addedPatchMap[patchI];
 515
 516             if (newPatchI != -1)
 517             {
 518                 const polyPatch& newPatch = mesh.boundaryMesh()[newPatchI];
 519                 const polyPatch& oldPatch =
 520                     fldToAdd.mesh().boundaryMesh()[patchI];
 521
 522                 if (!fld.boundaryField()(newPatchI))
 523                 {
 524                     // First occurrence of newPatchI. Map from existing
 525                     // patchField
 526
 527                     // From new patch faces to patch faces on added mesh.
 528                     labelList newToAdded
 529                     (
 530                         calcPatchMap
 531                         (
 532                             oldPatch.start(),
 533                             oldPatch.size(),
 534                             meshMap.addedFaceMap(),
 535                             newPatch,
 536                             0                       // unmapped values
 537                         )
 538                     );
 539
 540                     directFvPatchFieldMapper patchMapper(newToAdded);
 541
 542                     fld.boundaryField().set
 543                     (
 544                         newPatchI,
 545                         fvsPatchField<Type>::New
 546                         (
 547                             fldToAdd.boundaryField()[patchI],// added field
 548                             mesh.boundary()[newPatchI],      // new fvPatch
 549                             fld.dimensionedInternalField(),  // new int. field
 550                             patchMapper                      // mapper
 551                         )
 552                     );
 553                 }
 554                 else
 555                 {
 556                     // PatchField will have correct size already. Just slot in
 557                     // my elements.
 558
 559                     // From new patch faces to patch faces on added mesh. This
 560                     // time keep unmapped elements -1.
 561                     labelList newToAdded
 562                     (
 563                         calcPatchMap
 564                         (
 565                             oldPatch.start(),
 566                             oldPatch.size(),
 567                             meshMap.addedFaceMap(),
 568                             newPatch,
 569                             -1                      // unmapped values
 570                         )
 571                     );
 572
 573                     const fvsPatchField<Type>& addedFld =
 574                         fldToAdd.boundaryField()[patchI];
 575
 576                     fvsPatchField<Type>& newFld =
 577                         fld.boundaryField()[newPatchI];
 578
 579                     forAll(newFld, i)
 580                     {
 581                         label oldFaceI = newToAdded[i];
 582
 583                         if (oldFaceI >= 0 && oldFaceI < addedFld.size())
 584                         {
 585                             newFld[i] = addedFld[oldFaceI];
 586                         }
 587                     }
 588                 }
 589             }
 590         }
 591     }
 592 }
 593
 594
 595 template<class Type>
 596 void Foam::fvMeshAdder::MapSurfaceFields
 597 (
 598     const mapAddedPolyMesh& meshMap,
 599     const fvMesh& mesh,
 600     const fvMesh& meshToAdd
 601 )
 602 {
 603     typedef GeometricField<Type, fvsPatchField, surfaceMesh> fldType;
 604
 605     HashTable<const fldType*> fields
 606     (
 607         mesh.objectRegistry::lookupClass<fldType>()
 608     );
 609
 610     HashTable<const fldType*> fieldsToAdd
 611     (
 612         meshToAdd.objectRegistry::lookupClass<fldType>()
 613     );
 614
 615     // It is necessary to enforce that all old-time fields are stored
 616     // before the mapping is performed.  Otherwise, if the
 617     // old-time-level field is mapped before the field itself, sizes
 618     // will not match.
 619
 620     for
 621     (
 622         typename HashTable<const fldType*>::
 623             iterator fieldIter = fields.begin();
 624         fieldIter != fields.end();
 625         ++fieldIter
 626     )
 627     {
 628         const_cast<fldType*>(fieldIter())
 629             ->storeOldTimes();
 630     }
 631
 632
 633     for
 634     (
 635         typename HashTable<const fldType*>::
 636             iterator fieldIter = fields.begin();
 637         fieldIter != fields.end();
 638         ++fieldIter
 639     )
 640     {
 641         fldType& fld = const_cast<fldType&>(*fieldIter());
 642
 643         if (fieldsToAdd.found(fld.name()))
 644         {
 645             Pout<< "Mapping field " << fld.name() << endl;
 646
 647             const fldType& fldToAdd = *fieldsToAdd[fld.name()];
 648
 649             MapSurfaceField<Type>(meshMap, fld, fldToAdd);
 650         }
 651         else
 652         {
 653             WarningIn("fvMeshAdder::MapSurfaceFields")
 654                 << "Not mapping field " << fld.name()
 655                 << " since not present on mesh to add"
 656                 << endl;
 657         }
 658     }
 659 }
 660
 661
 662 // ************************************************************************* //