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