FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

findCellPointFaceTriangle.H

Go to the documentation of this file.
00001 /*
00002     find the triangle in which the position is,
00003     when the position is known to be on the face
00004 */
00005 
00006 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00007 
00008 namespace Foam
00009 {
00010 
00011 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00012 
00013 template<class Type>
00014 bool interpolationCellPointFace<Type>::findTriangle
00015 (
00016     const vector& position,
00017     const label nFace,
00018     label tetPointLabels[],
00019     scalar phi[]
00020 ) const
00021 {
00022     bool foundTriangle = false;
00023     vector tetPoints[3];
00024     const labelList& facePoints = this->pMeshFaces_[nFace];
00025     tetPoints[2] = this->pMeshFaceCentres_[nFace];
00026 
00027     label pointi = 0;
00028 
00029     while (pointi < facePoints.size() && !foundTriangle)
00030     {
00031         // The triangle is constructed from face center and one
00032         // face edge
00033         label nextPointLabel = (pointi + 1) % facePoints.size();
00034 
00035         tetPointLabels[0] = facePoints[pointi];
00036         tetPointLabels[1] = facePoints[nextPointLabel];
00037 
00038         tetPoints[0] = this->pMeshPoints_[tetPointLabels[0]];
00039         tetPoints[1] = this->pMeshPoints_[tetPointLabels[1]];
00040 
00041         vector fc = (tetPoints[0] + tetPoints[1] + tetPoints[2])/3.0;
00042 
00043         vector newPos = position + SMALL*(fc-position);
00044 
00045         // calculate triangle edge vectors and triangle face normal
00046         // the 'i':th edge is opposite node i
00047         vector edge[3], normal[3];
00048         for(label i=0; i<3; i++)
00049         {
00050             label ip0 = (i+1) % 3;
00051             label ipp = (i+2) % 3;
00052             edge[i] = tetPoints[ipp]-tetPoints[ip0];
00053         }
00054 
00055         vector triangleFaceNormal = edge[1] ^ edge[2];
00056         
00057         // calculate edge normal (pointing inwards)
00058         for(label i=0; i<3; i++)
00059         {
00060             normal[i] = triangleFaceNormal ^ edge[i];
00061             normal[i] /= mag(normal[i]) + VSMALL;
00062         }
00063 
00064         // check if position is inside triangle
00065         bool inside = true;
00066         for(label i=0; i<3; i++)
00067         {
00068             label ip = (i+1) % 3;
00069             inside = inside && (((newPos - tetPoints[ip]) & edge[i]) >= 0);
00070         }
00071 
00072         if (inside)
00073         {
00074             foundTriangle = true;
00075 
00076             // calculate phi-values
00077             for(label i=0; i<3; i++)
00078             {
00079                 label ip = (i+1) % 3;
00080                 scalar phiMax = max(VSMALL, normal[i] & edge[ip]);
00081                 scalar phiLength = (position-tetPoints[ip]) & normal[i];
00082                 phi[i] = phiLength/phiMax;
00083             }
00084         }
00085 
00086         pointi++;
00087     }
00088 
00089     return foundTriangle;
00090 }
00091 
00092 
00093 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00094 
00095 } // End namespace Foam
00096 
00097 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines