FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

Particle.C

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-2010 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
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your 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, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "Particle.H"
00027 #include <lagrangian/Cloud.H>
00028 #include <OpenFOAM/wedgePolyPatch.H>
00029 #include <OpenFOAM/symmetryPolyPatch.H>
00030 #include <OpenFOAM/cyclicPolyPatch.H>
00031 #include <OpenFOAM/processorPolyPatch.H>
00032 #include <OpenFOAM/wallPolyPatch.H>
00033 #include <OpenFOAM/transform.H>
00034 
00035 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00036 
00037 template<class ParticleType>
00038 void Foam::Particle<ParticleType>::findFaces
00039 (
00040     const vector& position,
00041     DynamicList<label>& faceList
00042 ) const
00043 {
00044     const polyMesh& mesh = cloud_.polyMesh_;
00045     const labelList& faces = mesh.cells()[celli_];
00046     const vector& C = mesh.cellCentres()[celli_];
00047 
00048     faceList.clear();
00049     forAll(faces, i)
00050     {
00051         label facei = faces[i];
00052         scalar lam = lambda(C, position, facei);
00053 
00054         if ((lam > 0) && (lam < 1.0))
00055         {
00056             faceList.append(facei);
00057         }
00058     }
00059 }
00060 
00061 
00062 template<class ParticleType>
00063 void Foam::Particle<ParticleType>::findFaces
00064 (
00065     const vector& position,
00066     const label celli,
00067     const scalar stepFraction,
00068     DynamicList<label>& faceList
00069 ) const
00070 {
00071     const polyMesh& mesh = cloud_.pMesh();
00072     const labelList& faces = mesh.cells()[celli];
00073     const vector& C = mesh.cellCentres()[celli];
00074 
00075     faceList.clear();
00076     forAll(faces, i)
00077     {
00078         label facei = faces[i];
00079         scalar lam = lambda(C, position, facei, stepFraction);
00080 
00081         if ((lam > 0) && (lam < 1.0))
00082         {
00083             faceList.append(facei);
00084         }
00085     }
00086 }
00087 
00088 
00089 template<class ParticleType>
00090 template<class TrackData>
00091 void Foam::Particle<ParticleType>::prepareForParallelTransfer
00092 (
00093     const label patchi,
00094     TrackData& td
00095 )
00096 {
00097     // Convert the face index to be local to the processor patch
00098     facei_ = patchFace(patchi, facei_);
00099 }
00100 
00101 
00102 template<class ParticleType>
00103 template<class TrackData>
00104 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
00105 (
00106     const label patchi,
00107     TrackData& td
00108 )
00109 {
00110     const processorPolyPatch& ppp =
00111         refCast<const processorPolyPatch>
00112         (cloud_.pMesh().boundaryMesh()[patchi]);
00113 
00114     celli_ = ppp.faceCells()[facei_];
00115 
00116     if (!ppp.parallel())
00117     {
00118         if (ppp.forwardT().size() == 1)
00119         {
00120             const tensor& T = ppp.forwardT()[0];
00121             transformPosition(T);
00122             static_cast<ParticleType&>(*this).transformProperties(T);
00123         }
00124         else
00125         {
00126             const tensor& T = ppp.forwardT()[facei_];
00127             transformPosition(T);
00128             static_cast<ParticleType&>(*this).transformProperties(T);
00129         }
00130     }
00131     else if (ppp.separated())
00132     {
00133         if (ppp.separation().size() == 1)
00134         {
00135             position_ -= ppp.separation()[0];
00136             static_cast<ParticleType&>(*this).transformProperties
00137             (
00138                 -ppp.separation()[0]
00139             );
00140         }
00141         else
00142         {
00143             position_ -= ppp.separation()[facei_];
00144             static_cast<ParticleType&>(*this).transformProperties
00145             (
00146                 -ppp.separation()[facei_]
00147             );
00148         }
00149     }
00150 
00151     // Reset the face index for the next tracking operation
00152     if (stepFraction_ > (1.0 - SMALL))
00153     {
00154         stepFraction_ = 1.0;
00155         facei_ = -1;
00156     }
00157     else
00158     {
00159         facei_ += ppp.start();
00160     }
00161 }
00162 
00163 
00164 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00165 
00166 template<class ParticleType>
00167 Foam::Particle<ParticleType>::Particle
00168 (
00169     const Cloud<ParticleType>& cloud,
00170     const vector& position,
00171     const label celli
00172 )
00173 :
00174     cloud_(cloud),
00175     position_(position),
00176     celli_(celli),
00177     facei_(-1),
00178     stepFraction_(0.0),
00179     origProc_(Pstream::myProcNo()),
00180     origId_(cloud_.getNewParticleID())
00181 {}
00182 
00183 
00184 template<class ParticleType>
00185 Foam::Particle<ParticleType>::Particle(const Particle<ParticleType>& p)
00186 :
00187     cloud_(p.cloud_),
00188     position_(p.position_),
00189     celli_(p.celli_),
00190     facei_(p.facei_),
00191     stepFraction_(p.stepFraction_),
00192     origProc_(p.origProc_),
00193     origId_(p.origId_)
00194 {}
00195 
00196 
00197 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00198 
00199 template<class ParticleType>
00200 template<class TrackData>
00201 Foam::label Foam::Particle<ParticleType>::track
00202 (
00203     const vector& endPosition,
00204     TrackData& td
00205 )
00206 {
00207     facei_ = -1;
00208 
00209     // Tracks to endPosition or stop on boundary
00210     while(!onBoundary() && stepFraction_ < 1.0 - SMALL)
00211     {
00212         stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
00213     }
00214 
00215     return facei_;
00216 }
00217 
00218 
00219 
00220 template<class ParticleType>
00221 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
00222 {
00223     int dummyTd;
00224     return track(endPosition, dummyTd);
00225 }
00226 
00227 template<class ParticleType>
00228 template<class TrackData>
00229 Foam::scalar Foam::Particle<ParticleType>::trackToFace
00230 (
00231     const vector& endPosition,
00232     TrackData& td
00233 )
00234 {
00235     const polyMesh& mesh = cloud_.polyMesh_;
00236 
00237     DynamicList<label>& faces = cloud_.labels_;
00238     findFaces(endPosition, faces);
00239 
00240     facei_ = -1;
00241     scalar trackFraction = 0.0;
00242 
00243     if (faces.empty())  // inside cell
00244     {
00245         trackFraction = 1.0;
00246         position_ = endPosition;
00247     }
00248     else // hit face
00249     {
00250         scalar lambdaMin = GREAT;
00251 
00252         if (faces.size() == 1)
00253         {
00254             lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
00255             facei_ = faces[0];
00256         }
00257         else
00258         {
00259             // If the particle has to cross more than one cell to reach the
00260             // endPosition, we check which way to go.
00261             // If one of the faces is a boundary face and the particle is
00262             // outside, we choose the boundary face.
00263             // The particle is outside if one of the lambda's is > 1 or < 0
00264             forAll(faces, i)
00265             {
00266                 scalar lam =
00267                     lambda(position_, endPosition, faces[i], stepFraction_);
00268 
00269                 if (lam < lambdaMin)
00270                 {
00271                     lambdaMin = lam;
00272                     facei_ = faces[i];
00273                 }
00274             }
00275         }
00276 
00277         bool internalFace = cloud_.internalFace(facei_);
00278 
00279         // For warped faces the particle can be 'outside' the cell.
00280         // This will yield a lambda larger than 1, or smaller than 0
00281         // For values < 0, the particle travels away from the cell
00282         // and we don't move the particle, only change cell.
00283         // For values larger than 1, we move the particle to endPosition only.
00284         if (lambdaMin > 0.0)
00285         {
00286             if (lambdaMin <= 1.0)
00287             {
00288                 trackFraction = lambdaMin;
00289                 position_ += trackFraction*(endPosition - position_);
00290             }
00291             else
00292             {
00293                 trackFraction = 1.0;
00294                 position_ = endPosition;
00295             }
00296         }
00297         else if (static_cast<ParticleType&>(*this).softImpact())
00298         {
00299             // Soft-sphere particles can travel outside the domain
00300             // but we don't use lambda since this the particle
00301             // is going away from face
00302             trackFraction = 1.0;
00303             position_ = endPosition;
00304         }
00305 
00306         // change cell
00307         if (internalFace) // Internal face
00308         {
00309             if (celli_ == mesh.faceOwner()[facei_])
00310             {
00311                 celli_ = mesh.faceNeighbour()[facei_];
00312             }
00313             else if (celli_ == mesh.faceNeighbour()[facei_])
00314             {
00315                 celli_ = mesh.faceOwner()[facei_];
00316             }
00317             else
00318             {
00319                 FatalErrorIn
00320                 (
00321                     "Particle::trackToFace(const vector&, TrackData&)"
00322                 )<< "addressing failure" << nl
00323                  << abort(FatalError);
00324             }
00325         }
00326         else
00327         {
00328             ParticleType& p = static_cast<ParticleType&>(*this);
00329 
00330             // Soft-sphere algorithm ignores the boundary
00331             if (p.softImpact())
00332             {
00333                 trackFraction = 1.0;
00334                 position_ = endPosition;
00335             }
00336 
00337             label patchi = patch(facei_);
00338             const polyPatch& patch = mesh.boundaryMesh()[patchi];
00339 
00340             if (!p.hitPatch(patch, td, patchi))
00341             {
00342                 if (isA<wedgePolyPatch>(patch))
00343                 {
00344                     p.hitWedgePatch
00345                     (
00346                         static_cast<const wedgePolyPatch&>(patch), td
00347                     );
00348                 }
00349                 else if (isA<symmetryPolyPatch>(patch))
00350                 {
00351                     p.hitSymmetryPatch
00352                     (
00353                         static_cast<const symmetryPolyPatch&>(patch), td
00354                     );
00355                 }
00356                 else if (isA<cyclicPolyPatch>(patch))
00357                 {
00358                     p.hitCyclicPatch
00359                     (
00360                         static_cast<const cyclicPolyPatch&>(patch), td
00361                     );
00362                 }
00363                 else if (isA<processorPolyPatch>(patch))
00364                 {
00365                     p.hitProcessorPatch
00366                     (
00367                         static_cast<const processorPolyPatch&>(patch), td
00368                     );
00369                 }
00370                 else if (isA<wallPolyPatch>(patch))
00371                 {
00372                     p.hitWallPatch
00373                     (
00374                         static_cast<const wallPolyPatch&>(patch), td
00375                     );
00376                 }
00377                 else
00378                 {
00379                     p.hitPatch(patch, td);
00380                 }
00381             }
00382         }
00383     }
00384 
00385     // If the trackFraction = 0 something went wrong.
00386     // Either the particle is flipping back and forth across a face perhaps
00387     // due to velocity interpolation errors or it is in a "hole" in the mesh
00388     // caused by face warpage.
00389     // In both cases resolve the positional ambiguity by moving the particle
00390     // slightly towards the cell-centre.
00391     if (trackFraction < SMALL)
00392     {
00393         position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
00394     }
00395 
00396     return trackFraction;
00397 }
00398 
00399 template<class ParticleType>
00400 Foam::scalar Foam::Particle<ParticleType>::trackToFace
00401 (
00402     const vector& endPosition
00403 )
00404 {
00405     int dummyTd;
00406     return trackToFace(endPosition, dummyTd);
00407 }
00408 
00409 template<class ParticleType>
00410 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
00411 {
00412     position_ = transform(T, position_);
00413 }
00414 
00415 
00416 template<class ParticleType>
00417 void Foam::Particle<ParticleType>::transformProperties(const tensor&)
00418 {}
00419 
00420 
00421 template<class ParticleType>
00422 void Foam::Particle<ParticleType>::transformProperties(const vector&)
00423 {}
00424 
00425 
00426 template<class ParticleType>
00427 template<class TrackData>
00428 bool Foam::Particle<ParticleType>::hitPatch
00429 (
00430     const polyPatch&,
00431     TrackData&,
00432     const label
00433 )
00434 {
00435     return false;
00436 }
00437 
00438 
00439 template<class ParticleType>
00440 template<class TrackData>
00441 void Foam::Particle<ParticleType>::hitWedgePatch
00442 (
00443     const wedgePolyPatch& wpp,
00444     TrackData&
00445 )
00446 {
00447     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
00448     nf /= mag(nf);
00449 
00450     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
00451 }
00452 
00453 
00454 template<class ParticleType>
00455 template<class TrackData>
00456 void Foam::Particle<ParticleType>::hitSymmetryPatch
00457 (
00458     const symmetryPolyPatch& spp,
00459     TrackData&
00460 )
00461 {
00462     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
00463     nf /= mag(nf);
00464 
00465     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
00466 }
00467 
00468 
00469 template<class ParticleType>
00470 template<class TrackData>
00471 void Foam::Particle<ParticleType>::hitCyclicPatch
00472 (
00473     const cyclicPolyPatch& cpp,
00474     TrackData&
00475 )
00476 {
00477     label patchFacei_ = cpp.whichFace(facei_);
00478 
00479     facei_ = cpp.transformGlobalFace(facei_);
00480 
00481     celli_ = cloud_.polyMesh_.faceOwner()[facei_];
00482 
00483     if (!cpp.parallel())
00484     {
00485         const tensor& T = cpp.transformT(patchFacei_);
00486 
00487         transformPosition(T);
00488         static_cast<ParticleType&>(*this).transformProperties(T);
00489     }
00490     else if (cpp.separated())
00491     {
00492         position_ += cpp.separation(patchFacei_);
00493         static_cast<ParticleType&>(*this).transformProperties
00494         (
00495             cpp.separation(patchFacei_)
00496         );
00497     }
00498 }
00499 
00500 
00501 template<class ParticleType>
00502 template<class TrackData>
00503 void Foam::Particle<ParticleType>::hitProcessorPatch
00504 (
00505     const processorPolyPatch& spp,
00506     TrackData& td
00507 )
00508 {}
00509 
00510 
00511 template<class ParticleType>
00512 template<class TrackData>
00513 void Foam::Particle<ParticleType>::hitWallPatch
00514 (
00515     const wallPolyPatch& spp,
00516     TrackData&
00517 )
00518 {}
00519 
00520 
00521 template<class ParticleType>
00522 template<class TrackData>
00523 void Foam::Particle<ParticleType>::hitPatch
00524 (
00525     const polyPatch& spp,
00526     TrackData&
00527 )
00528 {}
00529 
00530 
00531 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00532 
00533 #include "ParticleIO.C"
00534 
00535 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines