OpenFOAM logo
Open Source CFD Toolkit

particleI.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 
00027 \*---------------------------------------------------------------------------*/
00028 
00029 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00030 #include "polyMesh.H"
00031 
00032 namespace Foam
00033 {
00034 
00035 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
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         // move reference point for wall
00057         if (!cloud_.internalFace(facei))
00058         {
00059             const vector& C = mesh.cellCentres()[celli_];
00060             scalar CCf = mag((C - Cf) & Sf);
00061             // check if distance between cell centre and face centre
00062             // is larger than wallImpactDistance
00063             if (CCf > wallImpactDistance(Sf))
00064             {
00065                 Cf -= wallImpactDistance(Sf)*Sf;
00066             }
00067         }
00068         // for a moving mesh we need to reconstruct the old
00069         // Sf and Cf from oldPoints (they aren't stored)
00070         // NN.
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         // for layer addition Sf00 = vector::zero and we use Sf
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         // check if the face is rotating
00092         if (magSfDiff > SMALL)
00093         {
00094             vector Sf0 = Sf00 + fraction*(Sf - Sf00);
00095 
00096             // find center of rotation
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             // solve '(p - r0) & Sfp = 0', where
00105             // p = from + lambda*(to - from)
00106             // Sfp = Sf0 + lambda*(Sf - Sf0)
00107             // which results in the quadratic eq.
00108             // a*lambda^2 + b*lambda + c = 0
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                 // solve the second order polynomial
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                     // imaginary roots only
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                     // one root is around 0-1, while
00133                     // the other is very large in mag
00134                     if (mag(l1) < mag(l2))
00135                     {
00136                         return l1;
00137                     }
00138                     else
00139                     {
00140                         return l2;
00141                     }
00142                 }
00143             }
00144             else
00145             {
00146                 // when a==0, solve the first order polynomial
00147                 return (-c/b);
00148             }
00149         }
00150         else // no rotation
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             // check if trajectory is parallel to face
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         // mesh is static and fraction is not needed
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     // move reference point for wall
00196     if (!cloud_.internalFace(facei))
00197     {
00198         const vector& C = mesh.cellCentres()[celli_];
00199         scalar CCf = mag((C - Cf) & Sf);
00200         // check if distance between cell centre and face centre
00201         // is larger than wallImpactDistance
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     // check if trajectory is parallel to face
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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00249 
00250 //- Return reference to the particle cloud
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 //- Is the particle on a boundary face?
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 } // End namespace Foam
00299 
00300 // ************************************************************************* //
For further information go to www.openfoam.org