src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "volPointInterpolation.H"
  27 #include "volFields.H"
  28 #include "pointFields.H"
  29 #include "emptyFvPatch.H"
  30 #include "mapDistribute.H"
  31 #include "coupledPointPatchField.H"
  32 #include "valuePointPatchField.H"
  33 #include "transform.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 template<class Type, class CombineOp>
  43 void volPointInterpolation::syncUntransformedData
  44 (
  45     List<Type>& pointData,
  46     const CombineOp& cop
  47 ) const
  48 {
  49     // Transfer onto coupled patch
  50     const globalMeshData& gmd = mesh().globalData();
  51     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
  52     const labelList& meshPoints = cpp.meshPoints();
  53
  54     const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
  55     const labelListList& slaves = gmd.globalPointSlaves();
  56
  57     List<Type> elems(slavesMap.constructSize());
  58     forAll(meshPoints, i)
  59     {
  60         elems[i] = pointData[meshPoints[i]];
  61     }
  62
  63     // Pull slave data onto master. No need to update transformed slots.
  64     slavesMap.distribute(elems, false);
  65
  66     // Combine master data with slave data
  67     forAll(slaves, i)
  68     {
  69         Type& elem = elems[i];
  70
  71         const labelList& slavePoints = slaves[i];
  72
  73         // Combine master with untransformed slave data
  74         forAll(slavePoints, j)
  75         {
  76             cop(elem, elems[slavePoints[j]]);
  77         }
  78
  79         // Copy result back to slave slots
  80         forAll(slavePoints, j)
  81         {
  82             elems[slavePoints[j]] = elem;
  83         }
  84     }
  85
  86     // Push slave-slot data back to slaves
  87     slavesMap.reverseDistribute(elems.size(), elems, false);
  88
  89     // Extract back onto mesh
  90     forAll(meshPoints, i)
  91     {
  92         pointData[meshPoints[i]] = elems[i];
  93     }
  94 }
  95
  96
  97 template<class Type>
  98 void volPointInterpolation::pushUntransformedData
  99 (
 100     List<Type>& pointData
 101 ) const
 102 {
 103     // Transfer onto coupled patch
 104     const globalMeshData& gmd = mesh().globalData();
 105     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
 106     const labelList& meshPoints = cpp.meshPoints();
 107
 108     const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
 109     const labelListList& slaves = gmd.globalPointSlaves();
 110
 111     List<Type> elems(slavesMap.constructSize());
 112     forAll(meshPoints, i)
 113     {
 114         elems[i] = pointData[meshPoints[i]];
 115     }
 116
 117     // Combine master data with slave data
 118     forAll(slaves, i)
 119     {
 120         const labelList& slavePoints = slaves[i];
 121
 122         // Copy master data to slave slots
 123         forAll(slavePoints, j)
 124         {
 125             elems[slavePoints[j]] = elems[i];
 126         }
 127     }
 128
 129     // Push slave-slot data back to slaves
 130     slavesMap.reverseDistribute(elems.size(), elems, false);
 131
 132     // Extract back onto mesh
 133     forAll(meshPoints, i)
 134     {
 135         pointData[meshPoints[i]] = elems[i];
 136     }
 137 }
 138
 139
 140 template<class Type>
 141 void volPointInterpolation::addSeparated
 142 (
 143     GeometricField<Type, pointPatchField, pointMesh>& pf
 144 ) const
 145 {
 146     if (debug)
 147     {
 148         Pout<< "volPointInterpolation::addSeparated" << endl;
 149     }
 150
 151     forAll(pf.boundaryField(), patchI)
 152     {
 153         if (pf.boundaryField()[patchI].coupled())
 154         {
 155             refCast<coupledPointPatchField<Type> >
 156                 (pf.boundaryField()[patchI]).initSwapAddSeparated
 157                 (
 158                     Pstream::blocking,  //Pstream::nonBlocking,
 159                     pf.internalField()
 160                 );
 161         }
 162     }
 163
 164     // Block for any outstanding requests
 165     //Pstream::waitRequests();
 166
 167     forAll(pf.boundaryField(), patchI)
 168     {
 169         if (pf.boundaryField()[patchI].coupled())
 170         {
 171             refCast<coupledPointPatchField<Type> >
 172                 (pf.boundaryField()[patchI]).swapAddSeparated
 173                 (
 174                     Pstream::blocking,  //Pstream::nonBlocking,
 175                     pf.internalField()
 176                 );
 177         }
 178     }
 179 }
 180
 181
 182 template<class Type>
 183 void volPointInterpolation::interpolateInternalField
 184 (
 185     const GeometricField<Type, fvPatchField, volMesh>& vf,
 186     GeometricField<Type, pointPatchField, pointMesh>& pf
 187 ) const
 188 {
 189     if (debug)
 190     {
 191         Pout<< "volPointInterpolation::interpolateInternalField("
 192             << "const GeometricField<Type, fvPatchField, volMesh>&, "
 193             << "GeometricField<Type, pointPatchField, pointMesh>&) : "
 194             << "interpolating field from cells to points"
 195             << endl;
 196     }
 197
 198     const labelListList& pointCells = vf.mesh().pointCells();
 199
 200     // Multiply volField by weighting factor matrix to create pointField
 201     forAll(pointCells, pointi)
 202     {
 203         if (!isPatchPoint_[pointi])
 204         {
 205             const scalarList& pw = pointWeights_[pointi];
 206             const labelList& ppc = pointCells[pointi];
 207
 208             pf[pointi] = pTraits<Type>::zero;
 209
 210             forAll(ppc, pointCelli)
 211             {
 212                 pf[pointi] += pw[pointCelli]*vf[ppc[pointCelli]];
 213             }
 214         }
 215     }
 216 }
 217
 218
 219 template<class Type>
 220 tmp<Field<Type> > volPointInterpolation::flatBoundaryField
 221 (
 222     const GeometricField<Type, fvPatchField, volMesh>& vf
 223 ) const
 224 {
 225     const fvMesh& mesh = vf.mesh();
 226     const fvBoundaryMesh& bm = mesh.boundary();
 227
 228     tmp<Field<Type> > tboundaryVals
 229     (
 230         new Field<Type>(mesh.nFaces()-mesh.nInternalFaces())
 231     );
 232     Field<Type>& boundaryVals = tboundaryVals();
 233
 234     forAll(vf.boundaryField(), patchI)
 235     {
 236         label bFaceI = bm[patchI].patch().start() - mesh.nInternalFaces();
 237
 238         if (!isA<emptyFvPatch>(bm[patchI]) && !bm[patchI].coupled())
 239         {
 240             SubList<Type>
 241             (
 242                 boundaryVals,
 243                 vf.boundaryField()[patchI].size(),
 244                 bFaceI
 245             ).assign(vf.boundaryField()[patchI]);
 246         }
 247         else
 248         {
 249             const polyPatch& pp = bm[patchI].patch();
 250
 251             forAll(pp, i)
 252             {
 253                 boundaryVals[bFaceI++] = pTraits<Type>::zero;
 254             }
 255         }
 256     }
 257
 258     return tboundaryVals;
 259 }
 260
 261
 262 template<class Type>
 263 void volPointInterpolation::interpolateBoundaryField
 264 (
 265     const GeometricField<Type, fvPatchField, volMesh>& vf,
 266     GeometricField<Type, pointPatchField, pointMesh>& pf,
 267     const bool overrideFixedValue
 268 ) const
 269 {
 270     const primitivePatch& boundary = boundaryPtr_();
 271
 272     Field<Type>& pfi = pf.internalField();
 273
 274     // Get face data in flat list
 275     tmp<Field<Type> > tboundaryVals(flatBoundaryField(vf));
 276     const Field<Type>& boundaryVals = tboundaryVals();
 277
 278
 279     // Do points on 'normal' patches from the surrounding patch faces
 280     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 281
 282     forAll(boundary.meshPoints(), i)
 283     {
 284         label pointI = boundary.meshPoints()[i];
 285
 286         if (isPatchPoint_[pointI])
 287         {
 288             const labelList& pFaces = boundary.pointFaces()[i];
 289             const scalarList& pWeights = boundaryPointWeights_[i];
 290
 291             Type& val = pfi[pointI];
 292
 293             val = pTraits<Type>::zero;
 294             forAll(pFaces, j)
 295             {
 296                 if (boundaryIsPatchFace_[pFaces[j]])
 297                 {
 298                     val += pWeights[j]*boundaryVals[pFaces[j]];
 299                 }
 300             }
 301         }
 302     }
 303
 304     // Sum collocated contributions
 305     //mesh().globalData().syncPointData(pfi, plusEqOp<Type>());
 306     syncUntransformedData(pfi, plusEqOp<Type>());
 307
 308     // And add separated contributions
 309     addSeparated(pf);
 310
 311     // Push master data to slaves. It is possible (not sure how often) for
 312     // a coupled point to have its master on a different patch so
 313     // to make sure just push master data to slaves. Reuse the syncPointData
 314     // structure.
 315     //mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
 316     pushUntransformedData(pfi);
 317
 318
 319
 320     if (overrideFixedValue)
 321     {
 322         forAll(pf.boundaryField(), patchI)
 323         {
 324             pointPatchField<Type>& ppf = pf.boundaryField()[patchI];
 325
 326             if (isA<valuePointPatchField<Type> >(ppf))
 327             {
 328                 refCast<valuePointPatchField<Type> >(ppf) =
 329                     ppf.patchInternalField();
 330             }
 331         }
 332     }
 333
 334
 335     // Override constrained pointPatchField types with the constraint value.
 336     // This relys on only constrained pointPatchField implementing the evaluate
 337     // function
 338     pf.correctBoundaryConditions();
 339
 340     // Sync any dangling points
 341     //mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
 342     pushUntransformedData(pfi);
 343
 344     // Apply multiple constraints on edge/corner points
 345     applyCornerConstraints(pf);
 346 }
 347
 348
 349 template<class Type>
 350 void volPointInterpolation::applyCornerConstraints
 351 (
 352     GeometricField<Type, pointPatchField, pointMesh>& pf
 353 ) const
 354 {
 355     forAll(patchPatchPointConstraintPoints_, pointi)
 356     {
 357         pf[patchPatchPointConstraintPoints_[pointi]] = transform
 358         (
 359             patchPatchPointConstraintTensors_[pointi],
 360             pf[patchPatchPointConstraintPoints_[pointi]]
 361         );
 362     }
 363 }
 364
 365
 366 template<class Type>
 367 void volPointInterpolation::interpolate
 368 (
 369     const GeometricField<Type, fvPatchField, volMesh>& vf,
 370     GeometricField<Type, pointPatchField, pointMesh>& pf
 371 ) const
 372 {
 373     if (debug)
 374     {
 375         Pout<< "volPointInterpolation::interpolate("
 376             << "const GeometricField<Type, fvPatchField, volMesh>&, "
 377             << "GeometricField<Type, pointPatchField, pointMesh>&) : "
 378             << "interpolating field from cells to points"
 379             << endl;
 380     }
 381
 382     interpolateInternalField(vf, pf);
 383
 384     // Interpolate to the patches preserving fixed value BCs
 385     interpolateBoundaryField(vf, pf, false);
 386 }
 387
 388
 389 template<class Type>
 390 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 391 volPointInterpolation::interpolate
 392 (
 393     const GeometricField<Type, fvPatchField, volMesh>& vf,
 394     const wordList& patchFieldTypes
 395 ) const
 396 {
 397     const pointMesh& pm = pointMesh::New(vf.mesh());
 398
 399     // Construct tmp<pointField>
 400     tmp<GeometricField<Type, pointPatchField, pointMesh> > tpf
 401     (
 402         new GeometricField<Type, pointPatchField, pointMesh>
 403         (
 404             IOobject
 405             (
 406                 "volPointInterpolate(" + vf.name() + ')',
 407                 vf.instance(),
 408                 pm.thisDb()
 409             ),
 410             pm,
 411             vf.dimensions(),
 412             patchFieldTypes
 413         )
 414     );
 415
 416     interpolateInternalField(vf, tpf());
 417
 418     // Interpolate to the patches overriding fixed value BCs
 419     interpolateBoundaryField(vf, tpf(), true);
 420
 421     return tpf;
 422 }
 423
 424
 425 template<class Type>
 426 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 427 volPointInterpolation::interpolate
 428 (
 429     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf,
 430     const wordList& patchFieldTypes
 431 ) const
 432 {
 433     // Construct tmp<pointField>
 434     tmp<GeometricField<Type, pointPatchField, pointMesh> > tpf =
 435         interpolate(tvf(), patchFieldTypes);
 436     tvf.clear();
 437     return tpf;
 438 }
 439
 440
 441 template<class Type>
 442 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 443 volPointInterpolation::interpolate
 444 (
 445     const GeometricField<Type, fvPatchField, volMesh>& vf
 446 ) const
 447 {
 448     const pointMesh& pm = pointMesh::New(vf.mesh());
 449
 450     tmp<GeometricField<Type, pointPatchField, pointMesh> > tpf
 451     (
 452         new GeometricField<Type, pointPatchField, pointMesh>
 453         (
 454             IOobject
 455             (
 456                 "volPointInterpolate(" + vf.name() + ')',
 457                 vf.instance(),
 458                 pm.thisDb()
 459             ),
 460             pm,
 461             vf.dimensions()
 462         )
 463     );
 464
 465     interpolateInternalField(vf, tpf());
 466     interpolateBoundaryField(vf, tpf(), false);
 467
 468     return tpf;
 469 }
 470
 471
 472 template<class Type>
 473 tmp<GeometricField<Type, pointPatchField, pointMesh> >
 474 volPointInterpolation::interpolate
 475 (
 476     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf
 477 ) const
 478 {
 479     // Construct tmp<pointField>
 480     tmp<GeometricField<Type, pointPatchField, pointMesh> > tpf =
 481         interpolate(tvf());
 482     tvf.clear();
 483     return tpf;
 484 }
 485
 486
 487 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 488
 489 } // End namespace Foam
 490
 491 // ************************************************************************* //