src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 #include "fvMeshSubset.H"
  28 #include "emptyFvsPatchField.H"
  29 #include "emptyPointPatchField.H"
  30 #include "emptyFvPatchFields.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  38
  39 template<class Type>
  40 tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
  41 (
  42     const GeometricField<Type, fvPatchField, volMesh>& vf,
  43     const fvMesh& sMesh,
  44     const labelList& patchMap,
  45     const labelList& cellMap,
  46     const labelList& faceMap
  47 )
  48 {
  49     // Create and map the internal-field values
  50     Field<Type> internalField(vf.internalField(), cellMap);
  51
  52     // Create and map the patch field values
  53     PtrList<fvPatchField<Type> > patchFields(patchMap.size());
  54
  55     forAll (patchFields, patchI)
  56     {
  57         // Set the first one by hand as it corresponds to the
  58         // exposed internal faces.  Additional interpolation can be put here
  59         // as necessary.
  60         if (patchMap[patchI] == -1)
  61         {
  62             patchFields.set
  63             (
  64                 patchI,
  65                 new emptyFvPatchField<Type>
  66                 (
  67                     sMesh.boundary()[patchI],
  68                     DimensionedField<Type, volMesh>::null()
  69                 )
  70             );
  71         }
  72         else
  73         {
  74             // Construct addressing
  75             const fvPatch& subPatch = sMesh.boundary()[patchI];
  76             const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchI]];
  77             label baseStart = basePatch.patch().start();
  78             label baseSize = basePatch.size();
  79
  80             labelList directAddressing(subPatch.size());
  81
  82             forAll(directAddressing, i)
  83             {
  84                 label baseFaceI = faceMap[subPatch.patch().start()+i];
  85
  86                 if (baseFaceI >= baseStart && baseFaceI < baseStart+baseSize)
  87                 {
  88                     directAddressing[i] = baseFaceI-baseStart;
  89                 }
  90                 else
  91                 {
  92                     // Mapped from internal face. Do what? Map from element
  93                     // 0 for now.
  94                     directAddressing[i] = 0;
  95                 }
  96             }
  97
  98             patchFields.set
  99             (
 100                 patchI,
 101                 fvPatchField<Type>::New
 102                 (
 103                     vf.boundaryField()[patchMap[patchI]],
 104                     sMesh.boundary()[patchI],
 105                     DimensionedField<Type, volMesh>::null(),
 106                     patchFieldSubset(directAddressing)
 107                 )
 108             );
 109
 110             // What to do with exposed internal faces if put into this patch?
 111         }
 112     }
 113
 114
 115     // Create the complete field from the pieces
 116     tmp<GeometricField<Type, fvPatchField, volMesh> > tresF
 117     (
 118         new GeometricField<Type, fvPatchField, volMesh>
 119         (
 120             IOobject
 121             (
 122                 "subset"+vf.name(),
 123                 sMesh.time().timeName(),
 124                 sMesh,
 125                 IOobject::NO_READ,
 126                 IOobject::NO_WRITE
 127             ),
 128             sMesh,
 129             vf.dimensions(),
 130             internalField,
 131             patchFields
 132         )
 133     );
 134
 135     return tresF;
 136 }
 137
 138
 139 template<class Type>
 140 tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
 141 (
 142     const GeometricField<Type, fvPatchField, volMesh>& vf
 143 ) const
 144 {
 145     return interpolate
 146     (
 147         vf,
 148         subMesh(),
 149         patchMap(),
 150         cellMap(),
 151         faceMap()
 152     );
 153 }
 154
 155
 156 template<class Type>
 157 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
 158 (
 159     const GeometricField<Type, fvsPatchField, surfaceMesh>& vf,
 160     const fvMesh& sMesh,
 161     const labelList& patchMap,
 162     const labelList& faceMap
 163 )
 164 {
 165     // Create and map the internal-field values
 166     Field<Type> internalField
 167     (
 168         vf.internalField(),
 169         SubList<label>
 170         (
 171             faceMap,
 172             sMesh.nInternalFaces()
 173         )
 174     );
 175
 176     // Create and map the patch field values
 177     PtrList<fvsPatchField<Type> > patchFields(patchMap.size());
 178
 179     forAll (patchFields, patchI)
 180     {
 181         // Set the first one by hand as it corresponds to the
 182         // exposed internal faces.  Additional interpolation can be put here
 183         // as necessary.
 184         if (patchMap[patchI] == -1)
 185         {
 186             patchFields.set
 187             (
 188                 patchI,
 189                 new emptyFvsPatchField<Type>
 190                 (
 191                     sMesh.boundary()[patchI],
 192                     DimensionedField<Type, surfaceMesh>::null()
 193                 )
 194             );
 195         }
 196         else
 197         {
 198             // Construct addressing
 199             const fvPatch& subPatch = sMesh.boundary()[patchI];
 200             const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchI]];
 201             label baseStart = basePatch.patch().start();
 202             label baseSize = basePatch.size();
 203
 204             labelList directAddressing(subPatch.size());
 205
 206             forAll(directAddressing, i)
 207             {
 208                 label baseFaceI = faceMap[subPatch.patch().start()+i];
 209
 210                 if (baseFaceI >= baseStart && baseFaceI < baseStart+baseSize)
 211                 {
 212                     directAddressing[i] = baseFaceI-baseStart;
 213                 }
 214                 else
 215                 {
 216                     // Mapped from internal face. Do what? Map from element
 217                     // 0 for now.
 218                     directAddressing[i] = 0;
 219                 }
 220             }
 221
 222             patchFields.set
 223             (
 224                 patchI,
 225                 fvsPatchField<Type>::New
 226                 (
 227                     vf.boundaryField()[patchMap[patchI]],
 228                     sMesh.boundary()[patchI],
 229                     DimensionedField<Type, surfaceMesh>::null(),
 230                     patchFieldSubset(directAddressing)
 231                 )
 232             );
 233         }
 234     }
 235
 236
 237     // Map exposed internal faces. Note: Only nessecary if exposed faces added
 238     // into existing patch but since we don't know that at this point...
 239     forAll(patchFields, patchI)
 240     {
 241         fvsPatchField<Type>& pfld = patchFields[patchI];
 242
 243         label meshFaceI = pfld.patch().patch().start();
 244
 245         forAll(pfld, i)
 246         {
 247             label oldFaceI = faceMap[meshFaceI++];
 248
 249             if (oldFaceI < vf.internalField().size())
 250             {
 251                 pfld[i] = vf.internalField()[oldFaceI];
 252             }
 253         }
 254     }
 255
 256     // Create the complete field from the pieces
 257     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tresF
 258     (
 259         new GeometricField<Type, fvsPatchField, surfaceMesh>
 260         (
 261             IOobject
 262             (
 263                 "subset"+vf.name(),
 264                 sMesh.time().timeName(),
 265                 sMesh,
 266                 IOobject::NO_READ,
 267                 IOobject::NO_WRITE
 268             ),
 269             sMesh,
 270             vf.dimensions(),
 271             internalField,
 272             patchFields
 273         )
 274     );
 275
 276     return tresF;
 277 }
 278
 279
 280 template<class Type>
 281 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
 282 (
 283     const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
 284 ) const
 285 {
 286     return interpolate
 287     (
 288         sf,
 289         subMesh(),
 290         patchMap(),
 291         faceMap()
 292     );
 293 }
 294
 295
 296 template<class Type>
 297 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 298 fvMeshSubset::interpolate
 299 (
 300     const GeometricField<Type, pointPatchField, pointMesh>& vf,
 301     const pointMesh& sMesh,
 302     const labelList& patchMap,
 303     const labelList& pointMap
 304 )
 305 {
 306     // Create and map the internal-field values
 307     Field<Type> internalField(vf.internalField(), pointMap);
 308
 309     // Create and map the patch field values
 310     PtrList<pointPatchField<Type> > patchFields(patchMap.size());
 311
 312     forAll (patchFields, patchI)
 313     {
 314         // Set the first one by hand as it corresponds to the
 315         // exposed internal faces.  Additional interpolation can be put here
 316         // as necessary.
 317         if (patchMap[patchI] == -1)
 318         {
 319             patchFields.set
 320             (
 321                 patchI,
 322                 new emptyPointPatchField<Type>
 323                 (
 324                     sMesh.boundary()[patchI],
 325                     DimensionedField<Type, pointMesh>::null()
 326                 )
 327             );
 328         }
 329         else
 330         {
 331             // Construct addressing
 332             const pointPatch& basePatch =
 333                 vf.mesh().boundary()[patchMap[patchI]];
 334
 335             const labelList& meshPoints = basePatch.meshPoints();
 336
 337             // Make addressing from mesh to patch point
 338             Map<label> meshPointMap(2*meshPoints.size());
 339             forAll(meshPoints, localI)
 340             {
 341                 meshPointMap.insert(meshPoints[localI], localI);
 342             }
 343
 344             // Find which subpatch points originate from which patch point
 345             const pointPatch& subPatch = sMesh.boundary()[patchI];
 346             const labelList& subMeshPoints = subPatch.meshPoints();
 347
 348             // If mapped from outside patch use point 0 for lack of better.
 349             labelList directAddressing(subPatch.size(), 0);
 350
 351             forAll(subMeshPoints, localI)
 352             {
 353                 // Get mesh point on original mesh.
 354                 label meshPointI = pointMap[subMeshPoints[localI]];
 355
 356                 Map<label>::const_iterator iter = meshPointMap.find(meshPointI);
 357
 358                 if (iter != meshPointMap.end())
 359                 {
 360                     directAddressing[localI] = iter();
 361                 }
 362             }
 363
 364             patchFields.set
 365             (
 366                 patchI,
 367                 pointPatchField<Type>::New
 368                 (
 369                     vf.boundaryField()[patchMap[patchI]],
 370                     subPatch,
 371                     DimensionedField<Type, pointMesh>::null(),
 372                     pointPatchFieldSubset(directAddressing)
 373                 )
 374             );
 375         }
 376     }
 377
 378     // Create the complete field from the pieces
 379     tmp<GeometricField<Type, pointPatchField, pointMesh> > tresF
 380     (
 381         new GeometricField<Type, pointPatchField, pointMesh>
 382         (
 383             IOobject
 384             (
 385                 "subset"+vf.name(),
 386                 vf.time().timeName(),
 387                 sMesh.thisDb(),
 388                 IOobject::NO_READ,
 389                 IOobject::NO_WRITE
 390             ),
 391             sMesh,
 392             vf.dimensions(),
 393             internalField,
 394             patchFields
 395         )
 396     );
 397
 398     return tresF;
 399 }
 400
 401
 402 template<class Type>
 403 tmp<GeometricField<Type, pointPatchField, pointMesh> > fvMeshSubset::interpolate
 404 (
 405     const GeometricField<Type, pointPatchField, pointMesh>& sf
 406 ) const
 407 {
 408     return interpolate
 409     (
 410         sf,
 411         pointMesh::New(subMesh()),     // subsetted point mesh
 412         patchMap(),
 413         pointMap()
 414     );
 415 }
 416
 417
 418 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 419
 420 } // End namespace Foam
 421
 422 // ************************************************************************* //