00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 namespace Foam
00037 {
00038
00039
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
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
00102
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 }
00149
00150