diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C index e8efde1166..d53a904f3b 100644 --- a/src/functionObjects/field/nearWallFields/findCellParticle.C +++ b/src/functionObjects/field/nearWallFields/findCellParticle.C @@ -105,7 +105,7 @@ bool Foam::findCellParticle::move while (td.keepParticle && !td.switchProcessor && stepFraction() < 1) { const scalar f = 1 - stepFraction(); - trackToFace(f*(end_ - start_), f, td); + trackToAndHitFace(f*(end_ - start_), f, td); } if (stepFraction() == 1 || !td.keepParticle) diff --git a/src/functionObjects/field/streamLine/streamLineParticle.C b/src/functionObjects/field/streamLine/streamLineParticle.C index c587a7279c..58dd0b4030 100644 --- a/src/functionObjects/field/streamLine/streamLineParticle.C +++ b/src/functionObjects/field/streamLine/streamLineParticle.C @@ -198,7 +198,7 @@ bool Foam::streamLineParticle::move dt = maxDt; } - trackToFace(dt*U, 0, td); + trackToAndHitFace(dt*U, 0, td); if ( diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C index 2265ca5f84..7df06d3b98 100644 --- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C +++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C @@ -59,7 +59,7 @@ bool Foam::DSMCParcel::move(TrackData& td, const scalar trackTime) meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking); const scalar f = 1 - p.stepFraction(); - p.trackToFace(f*trackTime*Utracking, f, td); + p.trackToAndHitFace(f*trackTime*Utracking, f, td); if (p.onBoundaryFace() && td.keepParticle) { diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C index 2cef6b8f8b..e0d2d8df9a 100644 --- a/src/lagrangian/basic/particle/particle.C +++ b/src/lagrangian/basic/particle/particle.C @@ -344,6 +344,7 @@ void Foam::particle::changeFace(const label tetTriI) { const label newFaceI = mesh_.cells()[celli_][cellFaceI]; const class face& newFace = mesh_.faces()[newFaceI]; + const label newOwner = mesh_.faceOwner()[newFaceI]; // Exclude the current face if (tetFacei_ == newFaceI) @@ -351,21 +352,21 @@ void Foam::particle::changeFace(const label tetTriI) continue; } - // Loop over the edges, looking for the shared one - label edgeComp = 0; + // Loop over the edges, looking for the shared one. Note that we have to + // match the direction of the edge as well as the end points in order to + // avoid false positives when dealing with coincident ACMI faces. + const label edgeComp = newOwner == celli_ ? -1 : +1; label edgeI = 0; - for (; edgeI < newFace.size(); ++ edgeI) - { - edgeComp = edge::compare(sharedEdge, newFace.faceEdge(edgeI)); - - if (edgeComp != 0) - { - break; - } - } + for + ( + ; + edgeI < newFace.size() + && edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp; + ++ edgeI + ); // If the face does not contain the edge, then move on to the next face - if (edgeComp == 0) + if (edgeI >= newFace.size()) { continue; } @@ -643,6 +644,17 @@ Foam::scalar Foam::particle::trackToFace { // The track has hit a face, so set the current face and return facei_ = tetFacei_; + + /* + if (!mesh_.isInternalFace(facei_)) + { + const label patchi = mesh_.boundaryMesh().whichPatch(facei_); + const polyPatch& patch = mesh_.boundaryMesh()[patchi]; + Info<< "Tet face is on " << patch.type() << " patch #" << patchi + << endl << endl; + } + */ + return f; } else @@ -667,8 +679,8 @@ Foam::scalar Foam::particle::trackToStationaryTri if (debug) { - Info<< "Particle " << origId() << endl - << "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl << "Tracking from " << x0 + << " along " << x1 << " to " << x0 + x1 << endl; } // Get the tet geometry @@ -778,8 +790,8 @@ Foam::scalar Foam::particle::trackToMovingTri if (debug) { - Info<< "Particle " << origId() << endl - << "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl << "Tracking from " << x0 + << " along " << x1 << " to " << x0 + x1 << endl; } // Get the tet geometry diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 0c21308758..3a6f8213ef 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -57,6 +57,8 @@ class particle; class polyPatch; class cyclicPolyPatch; +class cyclicAMIPolyPatch; +class cyclicACMIPolyPatch; class processorPolyPatch; class symmetryPlanePolyPatch; class symmetryPolyPatch; @@ -294,10 +296,6 @@ protected: // Patch interactions - //- Overridable function to handle the particle hitting a face - template - void hitFace(TrackData& td); - //- Overridable function to handle the particle hitting a // patch. Executed before other patch-hitting functions. // trackFraction is passed in to allow mesh motion to @@ -343,6 +341,16 @@ protected: const vector& direction ); + //- Overridable function to handle the particle hitting a + // cyclicACMIPatch + template + void hitCyclicACMIPatch + ( + const cyclicACMIPolyPatch&, + TrackData& td, + const vector& direction + ); + //- Overridable function to handle the particle hitting a // processorPatch template @@ -585,12 +593,21 @@ public: label& tetTriI ); - //- As non-templated particle::trackToFace, but with additional - // boundary handling. + //- Hit the current face. If the current face is internal than this + // crosses into the next cell. If it is a boundary face then this will + // interact the particle with the relevant patch. template - void trackToFace + void hitFace ( - const vector& displacement, + const vector& direction, + TrackData& td + ); + + //- Convenience function. Cobines trackToFace and hitFace + template + void trackToAndHitFace + ( + const vector& direction, const scalar fraction, TrackData& td ); diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index 27291a03e8..307dcafff4 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -27,6 +27,7 @@ License #include "cyclicPolyPatch.H" #include "cyclicAMIPolyPatch.H" +#include "cyclicACMIPolyPatch.H" #include "processorPolyPatch.H" #include "symmetryPlanePolyPatch.H" #include "symmetryPolyPatch.H" @@ -186,38 +187,30 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr) template -void Foam::particle::trackToFace +void Foam::particle::hitFace ( - const vector& displacement, - const scalar fraction, + const vector& direction, TrackData& td ) { - // Track - trackToFace(displacement, fraction); - - // If the track is complete, return if (!onFace()) { return; } - - // Hit face/patch processing - typedef typename TrackData::cloudType::particleType particleType; - particleType& p = static_cast(*this); - p.hitFace(td); - if (onInternalFace()) + else if (onInternalFace()) { changeCell(); } - else + else if (onBoundaryFace()) { - label origFacei = facei_; - label patchi = mesh_.boundaryMesh().whichPatch(facei_); + typename TrackData::cloudType::particleType& p = + static_cast(*this); // No action is taken for tetPti_ for tetFacei_ here. These are handled // by the patch interaction call or later during processor transfer. + const label origFacei = facei_; + label patchi = mesh_.boundaryMesh().whichPatch(facei_); const tetIndices faceHitTetIs(celli_, tetFacei_, tetPti_); if @@ -268,13 +261,22 @@ void Foam::particle::trackToFace static_cast(patch), td ); } + else if (isA(patch)) + { + p.hitCyclicACMIPatch + ( + static_cast(patch), + td, + direction + ); + } else if (isA(patch)) { p.hitCyclicAMIPatch ( static_cast(patch), td, - displacement + direction ); } else if (isA(patch)) @@ -301,8 +303,16 @@ void Foam::particle::trackToFace template -void Foam::particle::hitFace(TrackData&) -{} +void Foam::particle::trackToAndHitFace +( + const vector& direction, + const scalar fraction, + TrackData& td +) +{ + trackToFace(direction, fraction); + hitFace(direction, td); +} template @@ -478,6 +488,46 @@ void Foam::particle::hitCyclicAMIPatch } +template +void Foam::particle::hitCyclicACMIPatch +( + const cyclicACMIPolyPatch& cpp, + TrackData& td, + const vector& direction +) +{ + const label localFacei = cpp.whichFace(facei_); + + // If the mask is within the patch tolerance at either end, then we can + // assume an interaction with the appropriate part of the ACMI pair. + const scalar mask = cpp.mask()[localFacei]; + bool couple = mask >= 1 - cpp.tolerance(); + bool nonOverlap = mask <= cpp.tolerance(); + + // If the mask is an intermediate value, then we search for a location on + // the other side of the AMI. If we can't find a location, then we assume + // that we have hit the non-overlap patch. + if (!couple && !nonOverlap) + { + vector pos = position(); + couple = cpp.pointFace(localFacei, direction, pos) >= 0; + nonOverlap = !couple; + } + + if (couple) + { + hitCyclicAMIPatch(cpp, td, direction); + } + else + { + // Move to the face associated with the non-overlap patch and redo the + // face interaction. + tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei; + hitFace(direction, td); + } +} + + template void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&) {} diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 35f2021187..037ada186c 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -292,8 +292,8 @@ bool Foam::KinematicParcel::move f = min(f, maxCo*l/max(SMALL*l, mag(s))); if (p.active()) { - // Track to the next face - p.trackToFace(f*s, f, td); + // Track to and hit the next face + p.trackToAndHitFace(f*s, f, td); } else { @@ -305,7 +305,6 @@ bool Foam::KinematicParcel::move p.stepFraction() += f; } - const scalar dt = (p.stepFraction() - sfrac)*trackTime; // Avoid problems with extremely small timesteps @@ -332,6 +331,11 @@ bool Foam::KinematicParcel::move p.age() += dt; + if (p.onFace()) + { + td.cloud().functions().postFace(p, p.face(), td.keepParticle); + } + td.cloud().functions().postMove(p, celli, dt, start, td.keepParticle); } @@ -339,22 +343,6 @@ bool Foam::KinematicParcel::move } -template -template -void Foam::KinematicParcel::hitFace(TrackData& td) -{ - typename TrackData::cloudType::parcelType& p = - static_cast(*this); - - td.cloud().functions().postFace(p, p.face(), td.keepParticle); -} - - -template -void Foam::KinematicParcel::hitFace(int& td) -{} - - template template bool Foam::KinematicParcel::hitPatch diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H index a862c4f033..fbdd578e27 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H @@ -597,14 +597,6 @@ public: // Patch interactions - //- Overridable function to handle the particle hitting a face - // without trackData - void hitFace(int& td); - - //- Overridable function to handle the particle hitting a face - template - void hitFace(TrackData& td); - //- Overridable function to handle the particle hitting a patch // Executed before other patch-hitting functions template diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C index f565adbfac..ef11eaddaa 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C @@ -88,7 +88,7 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime) while (td.keepParticle && !td.switchProcessor && stepFraction() < 1) { const scalar f = 1 - stepFraction(); - trackToFace(f*trackTime*v_, f, td); + trackToAndHitFace(f*trackTime*v_, f, td); } } else if (td.part() == 2) diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C index 4b4b8034b9..80737165ec 100644 --- a/src/lagrangian/solidParticle/solidParticle.C +++ b/src/lagrangian/solidParticle/solidParticle.C @@ -58,7 +58,7 @@ bool Foam::solidParticle::move const scalar sfrac = stepFraction(); const scalar f = 1 - stepFraction(); - trackToFace(f*trackTime*U_, f, td); + trackToAndHitFace(f*trackTime*U_, f, td); const scalar dt = (stepFraction() - sfrac)*trackTime; diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C index 71f44569be..46c13d83a1 100644 --- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C +++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C @@ -139,7 +139,8 @@ bool Foam::trackedParticle::move td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_); const scalar f = 1 - stepFraction(); - trackToFace(f*(end_ - start_), f, td); + const vector s = end_ - start_; + trackToAndHitFace(f*s, f, td); } } diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C index 429d1c084b..5ce5eaf85c 100644 --- a/src/sampling/sampledSet/face/faceOnlySet.C +++ b/src/sampling/sampledSet/face/faceOnlySet.C @@ -62,7 +62,7 @@ bool Foam::faceOnlySet::trackToBoundary { point oldPoint = trackPt; - singleParticle.trackToFace(end_ - start_, 0, trackData); + singleParticle.trackToAndHitFace(end_ - start_, 0, trackData); trackPt = singleParticle.position();