src/lagrangian/basic/Particle/ParticleI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "polyMesh.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 template<class ParticleType>
  37 inline scalar Particle<ParticleType>::lambda
  38 (
  39     const vector& from,
  40     const vector& to,
  41     const label facei,
  42     const scalar stepFraction
  43 ) const
  44 {
  45     const polyMesh& mesh = cloud_.polyMesh_;
  46     bool movingMesh = mesh.moving();
  47
  48     if (movingMesh)
  49     {
  50         vector Sf = mesh.faceAreas()[facei];
  51         Sf /= mag(Sf);
  52         vector Cf = mesh.faceCentres()[facei];
  53
  54         // move reference point for wall
  55         if (!cloud_.internalFace(facei))
  56         {
  57             const vector& C = mesh.cellCentres()[celli_];
  58             scalar CCf = mag((C - Cf) & Sf);
  59             // check if distance between cell centre and face centre
  60             // is larger than wallImpactDistance
  61             if
  62             (
  63                 CCf
  64               > static_cast<const ParticleType&>(*this).wallImpactDistance(Sf)
  65             )
  66             {
  67                 Cf -= static_cast<const ParticleType&>(*this)
  68                     .wallImpactDistance(Sf)*Sf;
  69             }
  70         }
  71
  72         // for a moving mesh we need to reconstruct the old
  73         // Sf and Cf from oldPoints (they aren't stored)
  74
  75         const vectorField& oldPoints = mesh.oldPoints();
  76
  77         vector Cf00 = mesh.faces()[facei].centre(oldPoints);
  78         vector Cf0 = Cf00 + stepFraction*(Cf - Cf00);
  79
  80         vector Sf00 = mesh.faces()[facei].normal(oldPoints);
  81
  82         // for layer addition Sf00 = vector::zero and we use Sf
  83         if (mag(Sf00) > SMALL)
  84         {
  85             Sf00 /= mag(Sf00);
  86         }
  87         else
  88         {
  89             Sf00 = Sf;
  90         }
  91
  92         scalar magSfDiff = mag(Sf - Sf00);
  93
  94         // check if the face is rotating
  95         if (magSfDiff > SMALL)
  96         {
  97             vector Sf0 = Sf00 + stepFraction*(Sf - Sf00);
  98
  99             // find center of rotation
 100             vector omega = Sf0 ^ Sf;
 101             scalar magOmega = mag(omega);
 102             omega /= magOmega + SMALL;
 103             vector n0 = omega ^ Sf0;
 104             scalar lam = ((Cf - Cf0) & Sf)/(n0 & Sf);
 105             vector r0 = Cf0 + lam*n0;
 106
 107             // solve '(p - r0) & Sfp = 0', where
 108             // p = from + lambda*(to - from)
 109             // Sfp = Sf0 + lambda*(Sf - Sf0)
 110             // which results in the quadratic eq.
 111             // a*lambda^2 + b*lambda + c = 0
 112             vector alpha = from - r0;
 113             vector beta = to - from;
 114             scalar a = beta & (Sf - Sf0);
 115             scalar b = (alpha & (Sf - Sf0)) + (beta & Sf0);
 116             scalar c = alpha & Sf0;
 117
 118             if (mag(a) > SMALL)
 119             {
 120                 // solve the second order polynomial
 121                 scalar ap = b/a;
 122                 scalar bp = c/a;
 123                 scalar cp = ap*ap - 4.0*bp;
 124
 125                 if (cp < 0.0)
 126                 {
 127                     // imaginary roots only
 128                     return GREAT;
 129                 }
 130                 else
 131                 {
 132                     scalar l1 = -0.5*(ap - ::sqrt(cp));
 133                     scalar l2 = -0.5*(ap + ::sqrt(cp));
 134
 135                     // one root is around 0-1, while
 136                     // the other is very large in mag
 137                     if (mag(l1) < mag(l2))
 138                     {
 139                         return l1;
 140                     }
 141                     else
 142                     {
 143                         return l2;
 144                     }
 145                 }
 146             }
 147             else
 148             {
 149                 // when a==0, solve the first order polynomial
 150                 return (-c/b);
 151             }
 152         }
 153         else // no rotation
 154         {
 155             vector alpha = from - Cf0;
 156             vector beta = to - from - (Cf - Cf0);
 157             scalar lambdaNominator = alpha & Sf;
 158             scalar lambdaDenominator = beta & Sf;
 159
 160             // check if trajectory is parallel to face
 161             if (mag(lambdaDenominator) < SMALL)
 162             {
 163                 if (lambdaDenominator < 0.0)
 164                 {
 165                     lambdaDenominator = -SMALL;
 166                 }
 167                 else
 168                 {
 169                     lambdaDenominator = SMALL;
 170                 }
 171             }
 172
 173             return (-lambdaNominator/lambdaDenominator);
 174         }
 175     }
 176     else
 177     {
 178         // mesh is static and stepFraction is not needed
 179         return lambda(from, to, facei);
 180     }
 181 }
 182
 183
 184 template<class ParticleType>
 185 inline scalar Particle<ParticleType>::lambda
 186 (
 187     const vector& from,
 188     const vector& to,
 189     const label facei
 190 ) const
 191 {
 192     const polyMesh& mesh = cloud_.polyMesh_;
 193
 194     vector Sf = mesh.faceAreas()[facei];
 195     Sf /= mag(Sf);
 196     vector Cf = mesh.faceCentres()[facei];
 197
 198     // move reference point for wall
 199     if (!cloud_.internalFace(facei))
 200     {
 201         const vector& C = mesh.cellCentres()[celli_];
 202         scalar CCf = mag((C - Cf) & Sf);
 203         // check if distance between cell centre and face centre
 204         // is larger than wallImpactDistance
 205         if
 206         (
 207             CCf
 208           > static_cast<const ParticleType&>(*this).wallImpactDistance(Sf)
 209         )
 210         {
 211             Cf -= static_cast<const ParticleType&>(*this)
 212                 .wallImpactDistance(Sf)*Sf;
 213         }
 214     }
 215
 216     scalar lambdaNominator = (Cf - from) & Sf;
 217     scalar lambdaDenominator = (to - from) & Sf;
 218
 219     // check if trajectory is parallel to face
 220     if (mag(lambdaDenominator) < SMALL)
 221     {
 222         if (lambdaDenominator < 0.0)
 223         {
 224             lambdaDenominator = -SMALL;
 225         }
 226         else
 227         {
 228             lambdaDenominator = SMALL;
 229         }
 230     }
 231
 232     return lambdaNominator/lambdaDenominator;
 233 }
 234
 235
 236 template<class ParticleType>
 237 inline bool Particle<ParticleType>::inCell() const
 238 {
 239     labelList faces = findFaces(position_);
 240
 241     return (!faces.size());
 242 }
 243
 244
 245 template<class ParticleType>
 246 inline bool Particle<ParticleType>::inCell
 247 (
 248     const vector& position,
 249     const label celli,
 250     const scalar stepFraction
 251 ) const
 252 {
 253     labelList faces = findFaces(position, celli, stepFraction);
 254
 255     return (!faces.size());
 256 }
 257
 258
 259 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 260
 261 template<class ParticleType>
 262 inline Particle<ParticleType>::trackData::trackData
 263 (
 264     Cloud<ParticleType>& cloud
 265 )
 266 :
 267     cloud_(cloud)
 268 {}
 269
 270 template<class ParticleType>
 271 inline Cloud<ParticleType>& Particle<ParticleType>::trackData::cloud()
 272 {
 273     return cloud_;
 274 }
 275
 276
 277 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 278
 279 template<class ParticleType>
 280 inline const Cloud<ParticleType>& Particle<ParticleType>::cloud() const
 281 {
 282     return cloud_;
 283 }
 284
 285
 286 template<class ParticleType>
 287 inline const vector& Particle<ParticleType>::position() const
 288 {
 289     return position_;
 290 }
 291
 292
 293 template<class ParticleType>
 294 inline vector& Particle<ParticleType>::position()
 295 {
 296     return position_;
 297 }
 298
 299
 300 template<class ParticleType>
 301 inline label Particle<ParticleType>::cell() const
 302 {
 303     return celli_;
 304 }
 305
 306 template<class ParticleType>
 307 inline label& Particle<ParticleType>::cell()
 308 {
 309     return celli_;
 310 }
 311
 312
 313 template<class ParticleType>
 314 inline label Particle<ParticleType>::face() const
 315 {
 316     return facei_;
 317 }
 318
 319
 320 template<class ParticleType>
 321 inline bool Particle<ParticleType>::onBoundary() const
 322 {
 323     return facei_ != -1 && facei_ >= cloud_.pMesh().nInternalFaces();
 324 }
 325
 326
 327 template<class ParticleType>
 328 inline scalar& Particle<ParticleType>::stepFraction()
 329 {
 330     return stepFraction_;
 331 }
 332
 333
 334 template<class ParticleType>
 335 inline scalar Particle<ParticleType>::stepFraction() const
 336 {
 337     return stepFraction_;
 338 }
 339
 340
 341 template<class ParticleType>
 342 inline bool Particle<ParticleType>::softImpact() const
 343 {
 344     return false;
 345 }
 346
 347
 348 template<class ParticleType>
 349 inline label Particle<ParticleType>::patch(const label facei) const
 350 {
 351     return cloud_.facePatch(facei);
 352 }
 353
 354
 355 template<class ParticleType>
 356 inline label Particle<ParticleType>::patchFace
 357 (
 358     const label patchi,
 359     const label facei
 360 ) const
 361 {
 362     return cloud_.patchFace(patchi, facei);
 363 }
 364
 365
 366 template<class ParticleType>
 367 inline scalar Particle<ParticleType>::wallImpactDistance(const vector&) const
 368 {
 369     return 0.0;
 370 }
 371
 372
 373 template<class ParticleType>
 374 inline label Particle<ParticleType>::faceInterpolation() const
 375 {
 376     return facei_;
 377 }
 378
 379
 380 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 381
 382 } // End namespace Foam
 383
 384 // ************************************************************************* //