src/lagrangian/basic/Particle/ParticleI.H

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