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 #include "polyMesh.H"
00031
00032 namespace Foam
00033 {
00034
00035
00036
00037 template<class particleType>
00038 inline scalar particle<particleType>::lambda
00039 (
00040 const vector& from,
00041 const vector& to,
00042 const label facei,
00043 const scalar fraction
00044 )
00045 {
00046 const polyMesh& mesh = cloud_.polyMesh_;
00047 bool movingMesh = mesh.moving();
00048
00049 if (movingMesh)
00050 {
00051
00052 vector Sf = mesh.faceAreas()[facei];
00053 Sf /= mag(Sf);
00054 vector Cf = mesh.faceCentres()[facei];
00055
00056
00057 if (!cloud_.internalFace(facei))
00058 {
00059 const vector& C = mesh.cellCentres()[celli_];
00060 scalar CCf = mag((C - Cf) & Sf);
00061
00062
00063 if (CCf > wallImpactDistance(Sf))
00064 {
00065 Cf -= wallImpactDistance(Sf)*Sf;
00066 }
00067 }
00068
00069
00070
00071
00072 const vectorField& oldPoints = mesh.oldPoints();
00073
00074 vector Cf00 = mesh.faces()[facei].centre(oldPoints);
00075 vector Cf0 = Cf00 + fraction*(Cf - Cf00);
00076
00077 vector Sf00 = mesh.faces()[facei].normal(oldPoints);
00078
00079
00080 if (mag(Sf00) > SMALL)
00081 {
00082 Sf00 /= mag(Sf00);
00083 }
00084 else
00085 {
00086 Sf00 = Sf;
00087 }
00088
00089 scalar magSfDiff = mag(Sf - Sf00);
00090
00091
00092 if (magSfDiff > SMALL)
00093 {
00094 vector Sf0 = Sf00 + fraction*(Sf - Sf00);
00095
00096
00097 vector omega = Sf0 ^ Sf;
00098 scalar magOmega = mag(omega);
00099 omega /= magOmega+SMALL;
00100 vector n0 = omega ^ Sf0;
00101 scalar lam = ((Cf - Cf0) & Sf)/(n0 & Sf);
00102 vector r0 = Cf0 + lam*n0;
00103
00104
00105
00106
00107
00108
00109 vector alpha = from - r0;
00110 vector beta = to - from;
00111 scalar a = beta & (Sf - Sf0);
00112 scalar b = (alpha & (Sf - Sf0)) + (beta & Sf0);
00113 scalar c = alpha & Sf0;
00114
00115 if (mag(a) > SMALL)
00116 {
00117
00118 scalar ap = b/a;
00119 scalar bp = c/a;
00120 scalar cp = ap*ap - 4.0*bp;
00121
00122 if (cp < 0.0)
00123 {
00124
00125 return GREAT;
00126 }
00127 else
00128 {
00129 scalar l1 = -0.5*(ap - ::sqrt(cp));
00130 scalar l2 = -0.5*(ap + ::sqrt(cp));
00131
00132
00133
00134 if (mag(l1) < mag(l2))
00135 {
00136 return l1;
00137 }
00138 else
00139 {
00140 return l2;
00141 }
00142 }
00143 }
00144 else
00145 {
00146
00147 return (-c/b);
00148 }
00149 }
00150 else
00151 {
00152 vector alpha = from - Cf0;
00153 vector beta = to - from - (Cf - Cf0);
00154 scalar lambdaNominator = alpha & Sf;
00155 scalar lambdaDenominator = beta & Sf;
00156
00157
00158 if (mag(lambdaDenominator) < SMALL)
00159 {
00160 if (lambdaDenominator < 0.0)
00161 {
00162 lambdaDenominator = -SMALL;
00163 }
00164 else
00165 {
00166 lambdaDenominator = SMALL;
00167 }
00168 }
00169
00170 return (-lambdaNominator/lambdaDenominator);
00171 }
00172 }
00173 else
00174 {
00175
00176 return lambda(from, to, facei);
00177 }
00178
00179 }
00180
00181 template<class particleType>
00182 inline scalar particle<particleType>::lambda
00183 (
00184 const vector& from,
00185 const vector& to,
00186 const label facei
00187 )
00188 {
00189 const polyMesh& mesh = cloud_.polyMesh_;
00190
00191 vector Sf = mesh.faceAreas()[facei];
00192 Sf /= mag(Sf);
00193 vector Cf = mesh.faceCentres()[facei];
00194
00195
00196 if (!cloud_.internalFace(facei))
00197 {
00198 const vector& C = mesh.cellCentres()[celli_];
00199 scalar CCf = mag((C - Cf) & Sf);
00200
00201
00202 if (CCf > wallImpactDistance(Sf))
00203 {
00204 Cf -= wallImpactDistance(Sf)*Sf;
00205 }
00206 }
00207
00208 scalar lambdaNominator = (Cf - from) & Sf;
00209 scalar lambdaDenominator = (to - from) & Sf;
00210
00211
00212 if (mag(lambdaDenominator) < SMALL)
00213 {
00214 if (lambdaDenominator < 0.0)
00215 {
00216 lambdaDenominator = -SMALL;
00217 }
00218 else
00219 {
00220 lambdaDenominator = SMALL;
00221 }
00222 }
00223
00224 return lambdaNominator/lambdaDenominator;
00225 }
00226
00227 template<class particleType>
00228 inline bool particle<particleType>::inCell()
00229 {
00230 labelList faces = findFaces(position_);
00231
00232 return (!faces.size());
00233 }
00234
00235 template<class particleType>
00236 inline bool particle<particleType>::inCell
00237 (
00238 const vector& position,
00239 const label celli,
00240 const scalar fraction
00241 )
00242 {
00243 labelList faces = findFaces(position, celli, fraction);
00244
00245 return (!faces.size());
00246 }
00247
00248
00249
00250
00251 template<class particleType>
00252 inline const Cloud<particleType>& particle<particleType>::cloud() const
00253 {
00254 return cloud_;
00255 }
00256
00257
00258 template<class particleType>
00259 inline const vector& particle<particleType>::position() const
00260 {
00261 return position_;
00262 }
00263
00264
00265 template<class particleType>
00266 inline label particle<particleType>::cell() const
00267 {
00268 return celli_;
00269 }
00270
00271
00272
00273 template<class particleType>
00274 inline bool particle<particleType>::onBoundary() const
00275 {
00276 return onBoundary_;
00277 }
00278
00279 template<class particleType>
00280 inline label particle<particleType>::patch(const label facei) const
00281 {
00282 return cloud_.facePatch(facei);
00283 }
00284
00285 template<class particleType>
00286 inline label particle<particleType>::patchFace
00287 (
00288 const label patchi,
00289 const label facei
00290 ) const
00291 {
00292 return cloud_.patchFace(patchi, facei);
00293 }
00294
00295
00296
00297
00298 }
00299
00300