src/lagrangian/basic/Particle/ParticleI.H

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