src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.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 "fvMeshSubset.H"
  28 #include "calculatedFvPatchFields.H"
  29 #include "calculatedFvsPatchFields.H"
  30 #include "emptyFvPatchFields.H"
  31 #include "emptyFvsPatchFields.H"
  32 #include "calculatedPointPatchFields.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  40
  41 template<class Type>
  42 tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::meshToMesh
  43 (
  44     const GeometricField<Type, fvPatchField, volMesh>& vf,
  45     const fvMesh& sMesh,
  46     const labelList& patchMap,
  47     const labelList& cellMap,
  48     const labelList& faceMap
  49 )
  50 {
  51     // Create and map the internal-field values
  52     Field<Type> internalField(vf.internalField(), cellMap);
  53
  54     // Create and map the patch field values
  55     PtrList<fvPatchField<Type> > patchFields(patchMap.size());
  56
  57     forAll (patchFields, patchI)
  58     {
  59         // Set the first one by hand as it corresponds to the
  60         // exposed internal faces.  Additional interpolation can be put here
  61         // as necessary.  HJ, date deleted
  62         if (patchMap[patchI] == -1)
  63         {
  64             patchFields.set
  65             (
  66                 patchI,
  67                 new emptyFvPatchField<Type>
  68                 (
  69                     sMesh.boundary()[patchI],
  70                     DimensionedField<Type, volMesh>::null()
  71                 )
  72             );
  73         }
  74         else
  75         {
  76             // Construct addressing
  77             const fvPatch& subPatch = sMesh.boundary()[patchI];
  78             const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchI]];
  79             label baseStart = basePatch.patch().start();
  80             label baseSize = basePatch.size();
  81
  82             labelList directAddressing(subPatch.size());
  83
  84             forAll(directAddressing, i)
  85             {
  86                 label baseFaceI = faceMap[subPatch.patch().start()+i];
  87
  88                 if (baseFaceI >= baseStart && baseFaceI < baseStart+baseSize)
  89                 {
  90                     directAddressing[i] = baseFaceI-baseStart;
  91                 }
  92                 else
  93                 {
  94                     // Mapped from internal face. Do what? Map from element
  95                     // 0 for now.
  96                     directAddressing[i] = 0;
  97                 }
  98             }
  99
 100             patchFields.set
 101             (
 102                 patchI,
 103                 fvPatchField<Type>::New
 104                 (
 105                     vf.boundaryField()[patchMap[patchI]],
 106                     sMesh.boundary()[patchI],
 107                     DimensionedField<Type, volMesh>::null(),
 108                     patchFieldSubset
 109                     (
 110                         vf.boundaryField()[patchMap[patchI]].size(),
 111                         directAddressing
 112                     )
 113                 )
 114             );
 115
 116             // What to do with exposed internal faces if put into this patch?
 117         }
 118     }
 119
 120
 121     // Create the complete field from the pieces
 122     tmp<GeometricField<Type, fvPatchField, volMesh> > tresF
 123     (
 124         new GeometricField<Type, fvPatchField, volMesh>
 125         (
 126             IOobject
 127             (
 128                 "subset"+vf.name(),
 129                 sMesh.time().timeName(),
 130                 sMesh,
 131                 IOobject::NO_READ,
 132                 IOobject::NO_WRITE
 133             ),
 134             sMesh,
 135             vf.dimensions(),
 136             internalField,
 137             patchFields
 138         )
 139     );
 140
 141     return tresF;
 142 }
 143
 144
 145 template<class Type>
 146 tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
 147 (
 148     const GeometricField<Type, fvPatchField, volMesh>& vf
 149 ) const
 150 {
 151     // Get reference to the subset mesh
 152     const fvMesh& sMesh = subMesh();
 153
 154     // Create and map the internal-field values
 155     Field<Type> internalField(vf.internalField(), cellMap());
 156
 157     // Create and map the patch field values
 158     const labelList& pm = patchMap();
 159
 160     // Create and map the patch field values
 161     PtrList<fvPatchField<Type> > patchFields(pm.size());
 162
 163     label internalFacesPatchIndex = -1;
 164
 165     forAll (patchFields, patchI)
 166     {
 167         // Set the first one by hand as it corresponds to the
 168         // exposed internal faces.  Additional interpolation can be put here
 169         // as necessary.  HJ, date deleted
 170         if (pm[patchI] == -1)
 171         {
 172             // Bug fix. Zeljko Tukovic, 10/Mar/2010
 173             internalFacesPatchIndex = patchI;
 174
 175             patchFields.set
 176             (
 177                 patchI,
 178                 new calculatedFvPatchField<Type>
 179                 (
 180                     sMesh.boundary()[patchI],
 181                     DimensionedField<Type, volMesh>::null()
 182                 )
 183             );
 184         }
 185         else
 186         {
 187             patchFields.set
 188             (
 189                 patchI,
 190                 fvPatchField<Type>::New
 191                 (
 192                     vf.boundaryField()[pm[patchI]],
 193                     sMesh.boundary()[patchI],
 194                     DimensionedField<Type, volMesh>::null(),
 195                     patchFieldSubset(*this, patchI)
 196                 )
 197             );
 198         }
 199     }
 200
 201     // Linear interpolation for last patch
 202     if (internalFacesPatchIndex > -1)
 203     {
 204         const Field<Type>& vfI = vf.internalField();
 205         const scalarField& w = baseMesh().weights().internalField();
 206         const labelList& own = baseMesh().faceOwner();
 207         const labelList& ngb = baseMesh().faceNeighbour();
 208
 209         Field<Type>& lastPatchField = patchFields[internalFacesPatchIndex];
 210
 211         label lastPatchStart =
 212             sMesh.boundaryMesh()[internalFacesPatchIndex].start();
 213
 214         const labelList& fm = faceMap();
 215
 216         forAll(lastPatchField, faceI)
 217         {
 218             lastPatchField[faceI] =
 219                 w[fm[lastPatchStart + faceI]]*
 220                 vfI[own[fm[lastPatchStart + faceI]]]
 221               + (1.0 - w[fm[lastPatchStart + faceI]])*
 222                 vfI[ngb[fm[lastPatchStart + faceI]]];
 223         }
 224     }
 225
 226     // Create the complete field from the pieces
 227     tmp<GeometricField<Type, fvPatchField, volMesh> > tresF
 228     (
 229         new GeometricField<Type, fvPatchField, volMesh>
 230         (
 231             IOobject
 232             (
 233                 "subset"+vf.name(),
 234                 sMesh.time().timeName(),
 235                 sMesh,
 236                 IOobject::NO_READ,
 237                 IOobject::NO_WRITE
 238             ),
 239             sMesh,
 240             vf.dimensions(),
 241             internalField,
 242             patchFields
 243         )
 244     );
 245
 246     return tresF;
 247 }
 248
 249
 250 template<class Type>
 251 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
 252 fvMeshSubset::meshToMesh
 253 (
 254     const GeometricField<Type, fvsPatchField, surfaceMesh>& vf,
 255     const fvMesh& sMesh,
 256     const labelList& patchMap,
 257     const labelList& faceMap
 258 )
 259 {
 260     // Create and map the internal-field values
 261     Field<Type> internalField
 262     (
 263         vf.internalField(),
 264         SubList<label>
 265         (
 266             faceMap,
 267             sMesh.nInternalFaces()
 268         )
 269     );
 270
 271     // Create and map the patch field values
 272     PtrList<fvsPatchField<Type> > patchFields(patchMap.size());
 273
 274     forAll (patchFields, patchI)
 275     {
 276         // Set the first one by hand as it corresponds to the
 277         // exposed internal faces.  Additional interpolation can be put here
 278         // as necessary.
 279         if (patchMap[patchI] == -1)
 280         {
 281             patchFields.set
 282             (
 283                 patchI,
 284                 new emptyFvsPatchField<Type>
 285                 (
 286                     sMesh.boundary()[patchI],
 287                     DimensionedField<Type, surfaceMesh>::null()
 288                 )
 289             );
 290         }
 291         else
 292         {
 293             // Construct addressing
 294             const fvPatch& subPatch = sMesh.boundary()[patchI];
 295             const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchI]];
 296             label baseStart = basePatch.patch().start();
 297             label baseSize = basePatch.size();
 298
 299             labelList directAddressing(subPatch.size());
 300
 301             forAll(directAddressing, i)
 302             {
 303                 label baseFaceI = faceMap[subPatch.patch().start()+i];
 304
 305                 if (baseFaceI >= baseStart && baseFaceI < baseStart+baseSize)
 306                 {
 307                     directAddressing[i] = baseFaceI-baseStart;
 308                 }
 309                 else
 310                 {
 311                     // Mapped from internal face. Do what? Map from element
 312                     // 0 for now.
 313                     directAddressing[i] = 0;
 314                 }
 315             }
 316
 317             patchFields.set
 318             (
 319                 patchI,
 320                 fvsPatchField<Type>::New
 321                 (
 322                     vf.boundaryField()[patchMap[patchI]],
 323                     sMesh.boundary()[patchI],
 324                     DimensionedField<Type, surfaceMesh>::null(),
 325                     patchFieldSubset
 326                     (
 327                         vf.boundaryField()[patchMap[patchI]].size(),
 328                         directAddressing
 329                     )
 330                 )
 331             );
 332         }
 333     }
 334
 335
 336     // Map exposed internal faces. Note: Only nessecary if exposed faces added
 337     // into existing patch but since we don't know that at this point...
 338     forAll(patchFields, patchI)
 339     {
 340         fvsPatchField<Type>& pfld = patchFields[patchI];
 341
 342         label meshFaceI = pfld.patch().patch().start();
 343
 344         forAll(pfld, i)
 345         {
 346             label oldFaceI = faceMap[meshFaceI++];
 347
 348             if (oldFaceI < vf.internalField().size())
 349             {
 350                 pfld[i] = vf.internalField()[oldFaceI];
 351             }
 352         }
 353     }
 354
 355     // Create the complete field from the pieces
 356     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tresF
 357     (
 358         new GeometricField<Type, fvsPatchField, surfaceMesh>
 359         (
 360             IOobject
 361             (
 362                 "subset"+vf.name(),
 363                 sMesh.time().timeName(),
 364                 sMesh,
 365                 IOobject::NO_READ,
 366                 IOobject::NO_WRITE
 367             ),
 368             sMesh,
 369             vf.dimensions(),
 370             internalField,
 371             patchFields
 372         )
 373     );
 374
 375     return tresF;
 376 }
 377
 378
 379 template<class Type>
 380 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
 381 fvMeshSubset::interpolate
 382 (
 383     const GeometricField<Type, fvsPatchField, surfaceMesh>& vf
 384 ) const
 385 {
 386     // Get reference to the subset mesh
 387     const fvMesh& sMesh = subMesh();
 388
 389     // Create and map the internal-field values
 390     Field<Type> internalField
 391     (
 392         vf.internalField(),
 393         SubList<label>
 394         (
 395             faceMap(),
 396             sMesh.nInternalFaces()
 397         )
 398     );
 399
 400     // Create and map the patch field values
 401     const labelList& pm = patchMap();
 402
 403     // Create and map the patch field values
 404     PtrList<fvsPatchField<Type> > patchFields(pm.size());
 405
 406     forAll (patchFields, patchI)
 407     {
 408         // Set the first one by hand as it corresponds to the
 409         // exposed internal faces.  Additional interpolation can be put here
 410         // as necessary.  HJ, date deleted
 411         if (pm[patchI] == -1)
 412         {
 413             patchFields.set
 414             (
 415                 patchI,
 416                 calculatedFvsPatchField<Type>
 417                 (
 418                     sMesh.boundary()[patchI],
 419                     DimensionedField<Type, surfaceMesh>::null()
 420                 )
 421             );
 422         }
 423         else
 424         {
 425             patchFields.set
 426             (
 427                 patchI,
 428                 fvsPatchField<Type>::New
 429                 (
 430                     vf.boundaryField()[pm[patchI]],
 431                     sMesh.boundary()[patchI],
 432                     DimensionedField<Type, surfaceMesh>::null(),
 433                     patchFieldSubset(*this, patchI)
 434                 )
 435             );
 436         }
 437     }
 438
 439
 440     const labelList& fm = faceMap();
 441
 442     // Map exposed internal faces. Note: Only nessecary if exposed faces added
 443     // into existing patch but since we don't know that at this point...
 444     forAll(patchFields, patchI)
 445     {
 446         fvsPatchField<Type>& pfld = patchFields[patchI];
 447
 448         label meshFaceI = pfld.patch().patch().start();
 449
 450         forAll(pfld, i)
 451         {
 452             label oldFaceI = fm[meshFaceI++];
 453
 454             if (oldFaceI < vf.internalField().size())
 455             {
 456                 pfld[i] = vf.internalField()[oldFaceI];
 457             }
 458         }
 459     }
 460
 461     // Create the complete field from the pieces
 462     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tresF
 463     (
 464         new GeometricField<Type, fvsPatchField, surfaceMesh>
 465         (
 466             IOobject
 467             (
 468                 "subset"+vf.name(),
 469                 sMesh.time().timeName(),
 470                 sMesh,
 471                 IOobject::NO_READ,
 472                 IOobject::NO_WRITE
 473             ),
 474             sMesh,
 475             vf.dimensions(),
 476             internalField,
 477             patchFields
 478         )
 479     );
 480
 481     return tresF;
 482 }
 483
 484
 485 template<class Type>
 486 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 487 fvMeshSubset::interpolate
 488 (
 489     const GeometricField<Type, pointPatchField, pointMesh>& vf
 490 ) const
 491 {
 492     // Get reference to the subset mesh
 493     const pointMesh& sMesh = subPointMesh();
 494
 495     // Create and map the internal-field values
 496     Field<Type> internalField(vf.internalField(), pointMap());
 497
 498     // Create and map the patch field values
 499     const labelList& pm = patchMap();
 500
 501     // Create and map the patch field values
 502     PtrList<pointPatchField<Type> > patchFields(pm.size());
 503
 504     forAll (patchFields, patchI)
 505     {
 506         // Set the first one by hand as it corresponds to the
 507         // exposed internal faces.  Additional interpolation can be put here
 508         // as necessary.  HJ, date deleted
 509         if (pm[patchI] == -1)
 510         {
 511             patchFields.set
 512             (
 513                 patchI,
 514                 new CalculatedPointPatchField
 515                 <
 516                     pointPatchField,
 517                     pointMesh,
 518                     pointPatch,
 519                     DummyMatrix,
 520                     Type
 521                 >
 522                 (
 523                     sMesh.boundary()[patchI],
 524                     DimensionedField<Type, pointMesh>::null()
 525                 )
 526             );
 527         }
 528         else
 529         {
 530             // Construct addressing
 531             const pointPatch& basePatch =
 532                 vf.mesh().boundary()[pm[patchI]];
 533
 534             const labelList& meshPoints = basePatch.meshPoints();
 535
 536             // Make addressing from mesh to patch point
 537             Map<label> meshPointMap(2*meshPoints.size());
 538             forAll(meshPoints, localI)
 539             {
 540                 meshPointMap.insert(meshPoints[localI], localI);
 541             }
 542
 543             // Find which subpatch points originate from which patch point
 544             const pointPatch& subPatch = sMesh.boundary()[patchI];
 545             const labelList& subMeshPoints = subPatch.meshPoints();
 546
 547             // If mapped from outside patch use point 0 for lack of better.
 548             labelList directAddressing(subPatch.size(), 0);
 549
 550             const labelList& ptMap = pointMap();
 551
 552             forAll(subMeshPoints, localI)
 553             {
 554                 // Get mesh point on original mesh.
 555                 label meshPointI = ptMap[subMeshPoints[localI]];
 556
 557                 Map<label>::const_iterator iter =
 558                     meshPointMap.find(meshPointI);
 559
 560                 if (iter != meshPointMap.end())
 561                 {
 562                     directAddressing[localI] = iter();
 563                 }
 564             }
 565
 566             patchFields.set
 567             (
 568                 patchI,
 569                 pointPatchField<Type>::New
 570                 (
 571                     vf.boundaryField()[pm[patchI]],
 572                     subPatch,
 573                     DimensionedField<Type, pointMesh>::null(),
 574                     pointPatchFieldSubset
 575                     (
 576                         directAddressing
 577                     )
 578                 )
 579             );
 580         }
 581     }
 582
 583     // Create the complete field from the pieces
 584     tmp<GeometricField<Type, pointPatchField, pointMesh> > tresF
 585     (
 586         new GeometricField<Type, pointPatchField, pointMesh>
 587         (
 588             IOobject
 589             (
 590                 "subset"+vf.name(),
 591                 vf.time().timeName(),
 592                 subMesh(),
 593                 IOobject::NO_READ,
 594                 IOobject::NO_WRITE
 595             ),
 596             sMesh,
 597             vf.dimensions(),
 598             internalField,
 599             patchFields
 600         )
 601     );
 602
 603     return tresF;
 604 }
 605
 606
 607 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 608
 609 } // End namespace Foam
 610
 611 // ************************************************************************* //