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