src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "interpolationCellPointFace.H"
  30 #include "volFields.H"
  31 #include "polyMesh.H"
  32 #include "volPointInterpolation.H"
  33 #include "linear.H"
  34 #include "findCellPointFaceTet.H"
  35 #include "findCellPointFaceTriangle.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
  43
  44 template<class Type>
  45 interpolationCellPointFace<Type>::interpolationCellPointFace
  46 (
  47     const volPointInterpolation& pInterp,
  48     const GeometricField<Type, fvPatchField, volMesh>& psi
  49 )
  50 :
  51     interpolation<Type>(psi),
  52     psip_(pInterp.interpolate(psi)),
  53     psis_(linearInterpolate(psi))
  54 {}
  55
  56
  57 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  58
  59 template<class Type>
  60 Type interpolationCellPointFace<Type>::interpolate
  61 (
  62     const vector& position,
  63     const label nCell,
  64     const label facei
  65 ) const
  66 {
  67     Type ts[4];
  68     vector tetPoints[4];
  69     scalar phi[4], phiCandidate[4];
  70     label tetLabelCandidate[2], tetPointLabels[2];
  71
  72     Type t = pTraits<Type>::zero;
  73
  74     // only use face information when the position is on a face
  75     if (facei < 0)
  76     {
  77         const vector& cellCentre = this->pMesh_.cellCentres()[nCell];
  78         const labelList& cellFaces = this->pMesh_.cells()[nCell];
  79
  80         vector projection = position - cellCentre;
  81         tetPoints[3] = cellCentre;
  82
  83         // *********************************************************************
  84         // project the cell-center through the point onto the face
  85         // and get the closest face, ...
  86         // *********************************************************************
  87
  88         bool foundTet = false;
  89         label closestFace = -1;
  90         scalar minDistance = GREAT;
  91
  92         forAll(cellFaces, facei)
  93         {
  94             label nFace = cellFaces[facei];
  95
  96             vector normal = this->pMeshFaceAreas_[nFace];
  97             normal /= mag(normal);
  98
  99             const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace];
 100
 101             scalar multiplierNumerator = (faceCentreTmp - cellCentre) & normal;
 102             scalar multiplierDenominator = projection & normal;
 103
 104             // if normal and projection are not orthogonal this could be the one...
 105             if (mag(multiplierDenominator) > SMALL)
 106             {
 107                 scalar multiplier = multiplierNumerator/multiplierDenominator;
 108                 vector iPoint = cellCentre + multiplier*projection;
 109                 scalar dist = mag(position - iPoint);
 110
 111                 if (dist < minDistance)
 112                 {
 113                     closestFace = nFace;
 114                     minDistance = dist;
 115                 }
 116             }
 117         }
 118
 119         // *********************************************************************
 120         // find the tetrahedron containing 'position'
 121         // from the cell center, face center and
 122         // two other points on the face
 123         // *********************************************************************
 124
 125         minDistance = GREAT;
 126         if (closestFace != -1)
 127         {
 128             label nFace = closestFace;
 129             foundTet = findTet
 130             (
 131                 position,
 132                 nFace,
 133                 tetPoints,
 134                 tetLabelCandidate,
 135                 tetPointLabels,
 136                 phi,
 137                 phiCandidate,
 138                 closestFace,
 139                 minDistance
 140             );
 141         }
 142
 143         if (!foundTet)
 144         {
 145             // check if the position is 'just' outside a tet
 146             if (minDistance < 1.0e-5)
 147             {
 148                 foundTet = true;
 149                 for (label i=0; i<4; i++)
 150                 {
 151                     phi[i] = phiCandidate[i];
 152                 }
 153                 tetPointLabels[0] = tetLabelCandidate[0];
 154                 tetPointLabels[1] = tetLabelCandidate[1];
 155             }
 156         }
 157
 158         // *********************************************************************
 159         // if the search failed check all the cell-faces
 160         // *********************************************************************
 161
 162         if (!foundTet)
 163         {
 164             minDistance = GREAT;
 165
 166             label facei = 0;
 167             while (facei < cellFaces.size() && !foundTet)
 168             {
 169                 label nFace = cellFaces[facei];
 170                 if (nFace < this->pMeshFaceAreas_.size())
 171                 {
 172                     foundTet = findTet
 173                     (
 174                         position,
 175                         nFace,
 176                         tetPoints,
 177                         tetLabelCandidate,
 178                         tetPointLabels,
 179                         phi,
 180                         phiCandidate,
 181                         closestFace,
 182                         minDistance
 183                     );
 184                 }
 185                 facei++;
 186             }
 187         }
 188
 189         if (!foundTet)
 190         {
 191             // check if the position is 'just' outside a tet
 192             // this time with a more tolerant limit
 193             if (minDistance < 1.0e-3)
 194             {
 195                 foundTet = true;
 196                 for (label i=0; i<4; i++)
 197                 {
 198                     phi[i] = phiCandidate[i];
 199                 }
 200                 tetPointLabels[0] = tetLabelCandidate[0];
 201                 tetPointLabels[1] = tetLabelCandidate[1];
 202             }
 203         }
 204
 205         // *********************************************************************
 206         // if the tet was found do the interpolation,
 207         // otherwise use the closest face value
 208         // *********************************************************************
 209
 210         if (foundTet)
 211         {
 212             for (label i=0; i<2; i++)
 213             {
 214                 ts[i] = psip_[tetPointLabels[i]];
 215             }
 216
 217             if (closestFace < psis_.size())
 218             {
 219                 ts[2] = psis_[closestFace];
 220             }
 221             else
 222             {
 223                 label patchi = this->pMesh_.boundaryMesh().whichPatch(closestFace);
 224
 225                 // If the boundary patch is not empty use the face value
 226                 // else use the cell value
 227                 if (this->psi_.boundaryField()[patchi].size())
 228                 {
 229                     ts[2] = this->psi_.boundaryField()[patchi]
 230                         [this->pMesh_.boundaryMesh()[patchi].whichFace(closestFace)];
 231                 }
 232                 else
 233                 {
 234                     ts[2] = this->psi_[nCell];
 235                 }
 236             }
 237
 238             ts[3] = this->psi_[nCell];
 239
 240             for (label n=0; n<4; n++)
 241             {
 242                 phi[n] = min(1.0, phi[n]);
 243                 phi[n] = max(0.0, phi[n]);
 244
 245                 t += phi[n]*ts[n];
 246             }
 247         }
 248         else
 249         {
 250             Info<< "interpolationCellPointFace<Type>::interpolate"
 251                 << "(const vector&, const label nCell) const : "
 252                 << "search failed; using closest cellFace value" << endl
 253                 << "cell number " << nCell << tab
 254                 << "position " << position << endl;
 255
 256             if (closestFace < psis_.size())
 257             {
 258                 t = psis_[closestFace];
 259             }
 260             else
 261             {
 262                 label patchi = this->pMesh_.boundaryMesh().whichPatch(closestFace);
 263
 264                 // If the boundary patch is not empty use the face value
 265                 // else use the cell value
 266                 if (this->psi_.boundaryField()[patchi].size())
 267                 {
 268                     t = this->psi_.boundaryField()[patchi]
 269                         [this->pMesh_.boundaryMesh()[patchi].whichFace(closestFace)];
 270                 }
 271                 else
 272                 {
 273                     t = this->psi_[nCell];
 274                 }
 275             }
 276         }
 277     }
 278     else
 279     {
 280         bool foundTriangle = findTriangle
 281         (
 282             position,
 283             facei,
 284             tetPointLabels,
 285             phi
 286         );
 287
 288         if (foundTriangle)
 289         {
 290             // add up the point values ...
 291             for (label i=0; i<2; i++)
 292             {
 293                 Type vel = psip_[tetPointLabels[i]];
 294                 t += phi[i]*vel;
 295             }
 296
 297             // ... and the face value
 298             if (facei < psis_.size())
 299             {
 300                 t += phi[2]*psis_[facei];
 301             }
 302             else
 303             {
 304                 label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
 305
 306                 // If the boundary patch is not empty use the face value
 307                 // else use the cell value
 308                 if (this->psi_.boundaryField()[patchi].size())
 309                 {
 310                     t += phi[2]*this->psi_.boundaryField()[patchi]
 311                         [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
 312                 }
 313                 else
 314                 {
 315                     t += phi[2]*this->psi_[nCell];
 316                 }
 317             }
 318         }
 319         else
 320         {
 321             // use face value only
 322             if (facei < psis_.size())
 323             {
 324                 t = psis_[facei];
 325             }
 326             else
 327             {
 328                 label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
 329
 330                 // If the boundary patch is not empty use the face value
 331                 // else use the cell value
 332                 if (this->psi_.boundaryField()[patchi].size())
 333                 {
 334                     t = this->psi_.boundaryField()[patchi]
 335                         [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
 336                 }
 337                 else
 338                 {
 339                     t = this->psi_[nCell];
 340                 }
 341             }
 342         }
 343     }
 344
 345     return t;
 346 }
 347
 348
 349 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 350
 351 } // End namespace Foam
 352
 353 // ************************************************************************* //