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