OpenFOAM logo
Open Source CFD Toolkit

findCellPointFaceTet.H

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software; you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by the
00013     Free Software Foundation; either version 2 of the License, or (at your
00014     option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM; if not, write to the Free Software Foundation,
00023     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
00024 
00025 Description
00026     find the tetrahedron in which the position is.
00027     while searching for the tet, store the tet
00028     closest to the position.
00029     This way, if position is not inside any tet,
00030     choose the closest one.
00031 
00032 \*---------------------------------------------------------------------------*/
00033 
00034 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00035 
00036 namespace Foam
00037 {
00038 
00039 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00040 
00041 template<class Type>
00042 bool interpolationCellPointFace<Type>::findTet
00043 (
00044     const vector& position,
00045     const label nFace,
00046     vector tetPoints[],
00047     label tetLabelCandidate[],
00048     label tetPointLabels[],
00049     scalar phi[],
00050     scalar phiCandidate[],
00051     label& closestFace,
00052     scalar& minDistance
00053 ) const
00054 {
00055     bool foundTet = false;
00056 
00057     const labelList& thisFacePoints = this->pMeshFaces_[nFace];
00058     tetPoints[2] = this->pMeshFaceCentres_[nFace];
00059 
00060     label pointi = 0;
00061 
00062     while (pointi < thisFacePoints.size() && !foundTet)
00063     {
00064         label nextPointLabel = (pointi + 1) % thisFacePoints.size();
00065 
00066         tetPointLabels[0] = thisFacePoints[pointi];
00067         tetPointLabels[1] = thisFacePoints[nextPointLabel];
00068 
00069         tetPoints[0] = this->pMeshPoints_[tetPointLabels[0]];
00070         tetPoints[1] = this->pMeshPoints_[tetPointLabels[1]];
00071 
00072         bool inside = true;
00073         scalar dist = 0.0;
00074 
00075         for (label n=0; n<4; n++)
00076         {
00077             label p1 = (n + 1) % 4;
00078             label p2 = (n + 2) % 4;
00079             label p3 = (n + 3) % 4;
00080 
00081             vector referencePoint, faceNormal;
00082             referencePoint = tetPoints[p1];
00083 
00084             faceNormal =
00085                 (tetPoints[p3] - tetPoints[p1])
00086               ^ (tetPoints[p2] - tetPoints[p1]);
00087 
00088             faceNormal /= mag(faceNormal);
00089 
00090             // correct normal to point into the tet
00091             vector v0 = tetPoints[n] - referencePoint;
00092             scalar correct = v0 & faceNormal;
00093             if (correct < 0)
00094             {
00095                 faceNormal = -faceNormal;
00096             }
00097 
00098             vector v1 = position - referencePoint + SMALL*faceNormal;
00099             scalar rightSide = v1 & faceNormal;
00100 
00101             // since normal is inwards, inside the tet
00102             // is defined as all dot-products being positive
00103             inside = inside && (rightSide >= 0);
00104 
00105             scalar phiLength = (position - referencePoint) & faceNormal;
00106 
00107             scalar maxLength =
00108                 max(VSMALL, (tetPoints[n] - referencePoint) & faceNormal);
00109 
00110             phi[n] = phiLength/maxLength;
00111 
00112             dist += phi[n];
00113         }
00114 
00115         if (!inside)
00116         {
00117             if (mag(dist - 1.0) < minDistance)
00118             {
00119                 minDistance = mag(dist - 1.0);
00120                 closestFace = nFace;
00121 
00122                 for (label i=0; i<4; i++)
00123                 {
00124                     phiCandidate[i] = phi[i];
00125                 }
00126 
00127                 tetLabelCandidate[0] = tetPointLabels[0];
00128                 tetLabelCandidate[1] = tetPointLabels[1];
00129             }
00130         }
00131 
00132         foundTet = inside;
00133 
00134         pointi++;
00135     }
00136 
00137     if (foundTet)
00138     {
00139         closestFace = nFace;
00140     }
00141 
00142     return foundTet;
00143 }
00144 
00145 
00146 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00147 
00148 } // End namespace Foam
00149 
00150 // ************************************************************************* //
For further information go to www.openfoam.org