diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C index 3df37ae160..2265ca5f84 100644 --- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C +++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C @@ -63,7 +63,7 @@ bool Foam::DSMCParcel::move(TrackData& td, const scalar trackTime) if (p.onBoundaryFace() && td.keepParticle) { - if (isA(pbMesh[p.patch(p.face())])) + if (isA(pbMesh[p.patch()])) { td.switchProcessor = true; } diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C index d206db410d..80d3683a00 100644 --- a/src/lagrangian/basic/particle/particle.C +++ b/src/lagrangian/basic/particle/particle.C @@ -40,7 +40,7 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::particle::tetGeometry +void Foam::particle::stationaryTetGeometry ( vector& centre, vector& base, @@ -59,23 +59,23 @@ void Foam::particle::tetGeometry } -Foam::barycentricTensor Foam::particle::tetTransform() const +Foam::barycentricTensor Foam::particle::stationaryTetTransform() const { vector centre, base, vertex1, vertex2; - tetGeometry(centre, base, vertex1, vertex2); + stationaryTetGeometry(centre, base, vertex1, vertex2); return barycentricTensor(centre, base, vertex1, vertex2); } -void Foam::particle::tetReverseTransform +void Foam::particle::stationaryTetReverseTransform ( vector& centre, scalar& detA, barycentricTensor& T ) const { - barycentricTensor A = tetTransform(); + barycentricTensor A = stationaryTetTransform(); const vector ab = A.b() - A.a(); const vector ac = A.c() - A.a(); @@ -117,27 +117,11 @@ void Foam::particle::movingTetGeometry const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces()); const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces()); - scalar f0 = stepFraction_, f1 = fraction; - // Old and new points and cell centres are not sub-cycled. If we are sub- // cycling, then we have to account for the timestep change here by // modifying the fractions that we take of the old and new geometry. - if (mesh_.time().subCycling()) - { - const TimeState& tsNew = mesh_.time(); - const TimeState& tsOld = mesh_.time().prevTimeState(); - - const scalar tFrac = - ( - (tsNew.value() - tsNew.deltaTValue()) - - (tsOld.value() - tsOld.deltaTValue()) - )/tsOld.deltaTValue(); - - const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue(); - - f0 = tFrac + dtFrac*f0; - f1 = dtFrac*f1; - } + const Pair s = stepFractionSpan(); + const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1]; centre[0] = ccOld + f0*(ccNew - ccOld); base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]); @@ -497,8 +481,8 @@ void Foam::particle::locate if (direction == nullptr) { - const polyPatch& p = mesh_.boundaryMesh()[patch(facei_)]; - direction = &p.faceNormals()[patchFace(patch(facei_), facei_)]; + const polyPatch& p = mesh_.boundaryMesh()[patch()]; + direction = &p.faceNormals()[p.whichFace(facei_)]; } const vector n = *direction/mag(*direction); @@ -646,14 +630,7 @@ Foam::scalar Foam::particle::trackToFace while (true) { - if (mesh_.moving()) - { - f *= trackToMovingTri(f*displacement, f*fraction, tetTriI); - } - else - { - f *= trackToStationaryTri(f*displacement, f*fraction, tetTriI); - } + f *= trackToTri(f*displacement, f*fraction, tetTriI); if (tetTriI == -1) { @@ -688,19 +665,20 @@ Foam::scalar Foam::particle::trackToStationaryTri if (debug) { - Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl + << "Tracking from " << x0 << " to " << x0 + x1 << endl; } // Get the tet geometry vector centre; scalar detA; barycentricTensor T; - tetReverseTransform(centre, detA, T); + stationaryTetReverseTransform(centre, detA, T); if (debug) { vector o, b, v1, v2; - tetGeometry(o, b, v1, v2); + stationaryTetGeometry(o, b, v1, v2); Info<< "Tet points o=" << o << ", b=" << b << ", v1=" << v1 << ", v2=" << v2 << endl << "Tet determinant = " << detA << endl @@ -717,7 +695,7 @@ Foam::scalar Foam::particle::trackToStationaryTri // Calculate the hit fraction label iH = -1; - scalar muH = detA == 0 ? VGREAT : mag(1/detA); + scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA; for (label i = 0; i < 4; ++ i) { if (std::isnormal(Tx1[i]) && Tx1[i] < 0) @@ -778,7 +756,7 @@ Foam::scalar Foam::particle::trackToStationaryTri // Set the proportion of the track that has been completed stepFraction_ += fraction*muH*detA; - return 1 - muH*detA; + return iH != -1 ? 1 - muH*detA : 0; } @@ -795,7 +773,8 @@ Foam::scalar Foam::particle::trackToMovingTri if (debug) { - Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl + << "Tracking from " << x0 << " to " << x0 + x1 << endl; } // Get the tet geometry @@ -823,7 +802,7 @@ Foam::scalar Foam::particle::trackToMovingTri // Calculate the hit fraction label iH = -1; - scalar muH = detA[0] == 0 ? VGREAT : mag(1/detA[0]); + scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0]; for (label i = 0; i < 4; ++ i) { const Roots<3> mu = hitEqn[i].roots(); @@ -899,7 +878,25 @@ Foam::scalar Foam::particle::trackToMovingTri << "End global coordinates = " << position() << endl; } - return 1 - muH*detA[0]; + return iH != -1 ? 1 - muH*detA[0] : 0; +} + + +Foam::scalar Foam::particle::trackToTri +( + const vector& displacement, + const scalar fraction, + label& tetTriI +) +{ + if (mesh_.moving()) + { + return trackToMovingTri(displacement, fraction, tetTriI); + } + else + { + return trackToStationaryTri(displacement, fraction, tetTriI); + } } @@ -926,6 +923,47 @@ void Foam::particle::constrainToMeshCentre() } +void Foam::particle::patchData(vector& n, vector& U) const +{ + if (!onBoundaryFace()) + { + FatalErrorInFunction + << "Patch data was requested for a particle that isn't on a patch" + << exit(FatalError); + } + + if (mesh_.moving()) + { + Pair centre, base, vertex1, vertex2; + movingTetGeometry(1, centre, base, vertex1, vertex2); + + n = triPointRef(base[0], vertex1[0], vertex2[0]).normal(); + n /= mag(n); + + // Interpolate the motion of the three face vertices to the current + // coordinates + U = + coordinates_.b()*base[1] + + coordinates_.c()*vertex1[1] + + coordinates_.d()*vertex2[1]; + + // The movingTetGeometry method gives the motion as a displacement + // across the time-step, so we divide by the time-step to get velocity + U /= mesh_.time().deltaTValue(); + } + else + { + vector centre, base, vertex1, vertex2; + stationaryTetGeometry(centre, base, vertex1, vertex2); + + n = triPointRef(base, vertex1, vertex2).normal(); + n /= mag(n); + + U = Zero; + } +} + + void Foam::particle::transformProperties(const tensor&) {} @@ -983,17 +1021,22 @@ void Foam::particle::correctAfterInteractionListReferral(const label celli) // position method is called. A referred particle should never be tracked, // so this approximate topology is good enough. By using the nearby cell we // minimise the error associated with the incorrect topology. + coordinates_ = barycentric(1, 0, 0, 0); if (mesh_.moving()) { - NotImplemented; + Pair centre; + FixedList detA; + FixedList T; + movingTetReverseTransform(0, centre, detA, T); + coordinates_ += (pos - centre[0]) & T[0]/detA[0]; } else { vector centre; scalar detA; barycentricTensor T; - tetReverseTransform(centre, detA, T); - coordinates_ = barycentric(1, 0, 0, 0) + ((pos - centre) & T/detA); + stationaryTetReverseTransform(centre, detA, T); + coordinates_ += (pos - centre) & T/detA; } } diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 2e3be525a7..7f6d486f79 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -175,7 +175,7 @@ private: // Tetrahedra functions //- Get the vertices of the current tet - void tetGeometry + void stationaryTetGeometry ( vector& centre, vector& base, @@ -188,7 +188,7 @@ private: // cartesian position in the global coordinate system. The // conversion is x = A & y, where x is the cartesian position, y is // the barycentric position and A is the transformation tensor. - barycentricTensor tetTransform() const; + barycentricTensor stationaryTetTransform() const; //- Get the reverse transform associated with the current tet. The // conversion is detA*y = (x - centre) & T. The variables x, y and @@ -197,7 +197,7 @@ private: // multiplied by its determinant, detA. This separation allows // the barycentric tracking algorithm to function on inverted or // degenerate tetrahedra. - void tetReverseTransform + void stationaryTetReverseTransform ( vector& centre, scalar& detA, @@ -454,10 +454,48 @@ public: //- Return current tet face particle is in inline label tetPt() const; + //- Return current face particle is on otherwise -1 + inline label face() const; + + //- Return the fraction of time-step completed + inline scalar stepFraction() const; + + //- Return the fraction of time-step completed + inline scalar& stepFraction(); + + //- Return the originating processor ID + inline label origProc() const; + + //- Return the originating processor ID + inline label& origProc(); + + //- Return the particle ID on the originating processor + inline label origId() const; + + //- Return the particle ID on the originating processor + inline label& origId(); + + + // Check + + //- Return the step fraction change within the overall time-step. + // Returns the start value and the change as a scalar pair. Always + // return Pair(0, 1), unless sub-cycling is in effect, in + // which case the values will reflect the span of the sub-cycle + // within the time-step. + inline Pair stepFractionSpan() const; + + //- Return the current fraction within the timestep. This differs + // from the stored step fraction due to sub-cycling. + inline scalar currentTimeFraction() const; + //- Return the indices of the current tet that the // particle occupies. inline tetIndices currentTetIndices() const; + //- Return the current tet transformation tensor + inline barycentricTensor currentTetTransform() const; + //- Return the normal of the tri on tetFacei_ for the // current tet. inline vector normal() const; @@ -467,15 +505,6 @@ public: // on oldPoints inline vector oldNormal() const; - //- Return current face particle is on otherwise -1 - inline label face() const; - - //- Return the particle current time - inline scalar currentTime() const; - - - // Check - //- Is the particle on a face? inline bool onFace() const; @@ -485,87 +514,80 @@ public: //- Is the particle on a boundary face? inline bool onBoundaryFace() const; - //- Which patch is particle on // <-- !!! - inline label patch(const label facei) const; - - //- Which face of this patch is this particle on // <-- !!! - inline label patchFace(const label patchi, const label facei) const; + //- Return the index of patch that the particle is on + inline label patch() const; //- Return current particle position inline vector position() const; - //- Return the fraction of time-step completed - inline scalar& stepFraction(); - //- Return the fraction of time-step completed - inline scalar stepFraction() const; + // Track - //- Return const access to the originating processor id - inline label origProc() const; + //- Track along the displacement for a given fraction of the overall + // step. End when the track is complete, or when a boundary is hit. + // On exit, stepFraction_ will have been incremented to the current + // position, and facei_ will be set to the index of the boundary + // face that was hit, or -1 if the track completed within a cell. + // The proportion of the displacement still to be completed is + // returned. + scalar track + ( + const vector& displacement, + const scalar fraction + ); - //- Return the originating processor id for manipulation - inline label& origProc(); + //- As particle::track, but also stops on internal faces. + scalar trackToFace + ( + const vector& displacement, + const scalar fraction + ); - //- Return const access to the particle id on originating processor - inline label origId() const; + //- As particle::trackToFace, but also stops on tet triangles. On + // exit, tetTriI is set to the index of the tet triangle that was + // hit, or -1 if the end position was reached. + scalar trackToTri + ( + const vector& displacement, + const scalar fraction, + label& tetTriI + ); - //- Return the particle id on originating processor for manipulation - inline label& origId(); + //- As particle::trackToTri, but for stationary meshes + scalar trackToStationaryTri + ( + const vector& displacement, + const scalar fraction, + label& tetTriI + ); + + //- As particle::trackToTri, but for moving meshes + scalar trackToMovingTri + ( + const vector& displacement, + const scalar fraction, + label& tetTriI + ); + + //- As non-templated particle::trackToFace, but with additional + // boundary handling. + template + void trackToFace + ( + const vector& displacement, + const scalar fraction, + TrackData& td + ); + + //- Set the constrained components of the particle position to the + // mesh centre. + void constrainToMeshCentre(); - // Track + // Patch data - //- Track along the displacement for a given fraction of the overall - // step. End when the track is complete, or when a boundary is hit. - // On exit, stepFraction_ will have been incremented to the current - // position, and facei_ will be set to the index of the boundary - // face that was hit, or -1 if the track completed within a cell. - // The proportion of the displacement still to be completed is - // returned. - scalar track - ( - const vector& displacement, - const scalar fraction - ); - - //- As particle::track, but also stops on internal faces. - scalar trackToFace - ( - const vector& displacement, - const scalar fraction - ); - - //- As particle::trackToFace, but also stops on tet triangles. On - // exit, tetTriI is set to the index of the tet triangle that was - // hit, or -1 if the end position was reached. - scalar trackToStationaryTri - ( - const vector& displacement, - const scalar fraction, - label& tetTriI - ); - - //- As particle::trackToTri, but for moving meshes - scalar trackToMovingTri - ( - const vector& displacement, - const scalar fraction, - label& tetTriI - ); - - //- As non-templated particle::trackToFace, but with additional - // boundary handling. - template - void trackToFace - ( - const vector& displacement, - const scalar fraction, - TrackData& td - ); - - //- Set the constrained components of the particle position to the - // mesh centre. - void constrainToMeshCentre(); + //- Get the normal and velocity of the current patch location + void patchData(vector& n, vector& U) const; // Transformations diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H index fce5ded593..7249fcc275 100644 --- a/src/lagrangian/basic/particle/particleI.H +++ b/src/lagrangian/basic/particle/particleI.H @@ -72,68 +72,19 @@ inline Foam::label Foam::particle::tetPt() const } -inline Foam::tetIndices Foam::particle::currentTetIndices() const -{ - return tetIndices(celli_, tetFacei_, tetPti_); -} - - -inline Foam::vector Foam::particle::normal() const -{ - return currentTetIndices().faceTri(mesh_).normal(); -} - - -inline Foam::vector Foam::particle::oldNormal() const -{ - return currentTetIndices().oldFaceTri(mesh_).normal(); -} - - inline Foam::label Foam::particle::face() const { return facei_; } -inline bool Foam::particle::onFace() const -{ - return facei_ >= 0; -} - - -inline bool Foam::particle::onInternalFace() const -{ - return onFace() && mesh_.isInternalFace(facei_); -} - - -inline bool Foam::particle::onBoundaryFace() const -{ - return onFace() && !mesh_.isInternalFace(facei_); -} - - -inline Foam::vector Foam::particle::position() const -{ - if (mesh_.moving()) - { - return movingTetTransform(0)[0] & coordinates_; - } - else - { - return tetTransform() & coordinates_; - } -} - - -inline Foam::scalar& Foam::particle::stepFraction() +inline Foam::scalar Foam::particle::stepFraction() const { return stepFraction_; } -inline Foam::scalar Foam::particle::stepFraction() const +inline Foam::scalar& Foam::particle::stepFraction() { return stepFraction_; } @@ -163,25 +114,96 @@ inline Foam::label& Foam::particle::origId() } -inline Foam::scalar Foam::particle::currentTime() const +inline Foam::Pair Foam::particle::stepFractionSpan() const { - return mesh_.time().value() + stepFraction_*mesh_.time().deltaTValue(); + if (mesh_.time().subCycling()) + { + const TimeState& tsNew = mesh_.time(); + const TimeState& tsOld = mesh_.time().prevTimeState(); + + const scalar tFrac = + ( + (tsNew.value() - tsNew.deltaTValue()) + - (tsOld.value() - tsOld.deltaTValue()) + )/tsOld.deltaTValue(); + + const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue(); + + return Pair(tFrac, dtFrac); + } + else + { + return Pair(0, 1); + } } -inline Foam::label Foam::particle::patch(const label facei) const +inline Foam::scalar Foam::particle::currentTimeFraction() const { - return mesh_.boundaryMesh().whichPatch(facei); + const Pair s = stepFractionSpan(); + + return s[0] + stepFraction_*s[1]; } -inline Foam::label Foam::particle::patchFace -( - const label patchi, - const label facei -) const +inline Foam::tetIndices Foam::particle::currentTetIndices() const { - return mesh_.boundaryMesh()[patchi].whichFace(facei); + return tetIndices(celli_, tetFacei_, tetPti_); +} + + +inline Foam::barycentricTensor Foam::particle::currentTetTransform() const +{ + if (mesh_.moving()) + { + return movingTetTransform(0)[0]; + } + else + { + return stationaryTetTransform(); + } +} + + +inline Foam::vector Foam::particle::normal() const +{ + return currentTetIndices().faceTri(mesh_).normal(); +} + + +inline Foam::vector Foam::particle::oldNormal() const +{ + return currentTetIndices().oldFaceTri(mesh_).normal(); +} + + +inline bool Foam::particle::onFace() const +{ + return facei_ >= 0; +} + + +inline bool Foam::particle::onInternalFace() const +{ + return onFace() && mesh_.isInternalFace(facei_); +} + + +inline bool Foam::particle::onBoundaryFace() const +{ + return onFace() && !mesh_.isInternalFace(facei_); +} + + +inline Foam::label Foam::particle::patch() const +{ + return mesh_.boundaryMesh().whichPatch(facei_); +} + + +inline Foam::vector Foam::particle::position() const +{ + return currentTetTransform() & coordinates_; } diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C index 8b0291bc02..80ad21e2ed 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C @@ -703,143 +703,39 @@ void Foam::KinematicCloud::patchData ( const parcelType& p, const polyPatch& pp, - const scalar trackFraction, - const tetIndices& tetIs, vector& nw, vector& Up ) const { - label patchi = pp.index(); - label patchFacei = pp.whichFace(p.face()); + p.patchData(nw, Up); - vector n = tetIs.faceTri(mesh_).normal(); - n /= mag(n); - - vector U = U_.boundaryField()[patchi][patchFacei]; - - // Unless the face is rotating, the required normal is n; - nw = n; - - if (!mesh_.moving()) + // If this is a wall patch, then there may be a non-zero tangential velocity + // component; the lid velocity in a lid-driven cavity case, for example. We + // want the particle to interact with this velocity, so we look it up in the + // velocity field and use it to set the wall-tangential component. + if (isA(pp)) { - // Only wall patches may have a non-zero wall velocity from - // the velocity field when the mesh is not moving. + const label patchi = pp.index(); + const label patchFacei = pp.whichFace(p.face()); - if (isA(pp)) + // We only want to use the boundary condition value onlyif it is set + // by the boundary condition. If the boundary values are extrapolated + // (e.g., slip conditions) then they represent the motion of the fluid + // just inside the domain rather than that of the wall itself. + if (U_.boundaryField()[patchi].fixesValue()) { - Up = U; + const vector Uw1 = U_.boundaryField()[patchi][patchFacei]; + const vector& Uw0 = + U_.oldTime().boundaryField()[patchi][patchFacei]; + + const scalar f = p.currentTimeFraction(); + + const vector Uw = Uw0 + f*(Uw1 - Uw0); + + const tensor nnw = nw*nw; + + Up = (nnw & Up) + Uw - (nnw & Uw); } - else - { - Up = Zero; - } - } - else - { - vector U00 = U_.oldTime().boundaryField()[patchi][patchFacei]; - - vector n00 = tetIs.oldFaceTri(mesh_).normal(); - - // Difference in normal over timestep - vector dn = Zero; - - if (mag(n00) > SMALL) - { - // If the old normal is zero (for example in layer - // addition) then use the current normal, meaning that the - // motion can only be translational, and dn remains zero, - // otherwise, calculate dn: - - n00 /= mag(n00); - - dn = n - n00; - } - - // Total fraction through the timestep of the motion, - // including stepFraction before the current tracking step - // and the current trackFraction - // i.e. - // let s = stepFraction, t = trackFraction - // Motion of x in time: - // |-----------------|---------|---------| - // x00 x0 xi x - // - // where xi is the correct value of x at the required - // tracking instant. - // - // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00 - // - // i.e. the motion covered by previous tracking portions - // within this timestep, and - // - // xi = x0 + t*(x - x0) - // = t*x + (1 - t)*x0 - // = t*x + (1 - t)*(s*x + (1 - s)*x00) - // = (s + t - s*t)*x + (1 - (s + t - s*t))*x00 - // - // let m = (s + t - s*t) - // - // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00); - // - // In the same form as before. - - scalar m = - p.stepFraction() - + trackFraction - - (p.stepFraction()*trackFraction); - - // When the mesh is moving, the velocity field on wall patches - // will contain the velocity associated with the motion of the - // mesh, in which case it is interpolated in time using m. - // For other patches the face velocity will need to be - // reconstructed from the face centre motion. - - const vector& Cf = mesh_.faceCentres()[p.face()]; - - vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints()); - - if (isA(pp)) - { - Up = U00 + m*(U - U00); - } - else - { - Up = (Cf - Cf00)/mesh_.time().deltaTValue(); - } - - if (mag(dn) > SMALL) - { - // Rotational motion, nw requires interpolation and a - // rotational velocity around face centre correction to Up - // is required. - - nw = n00 + m*dn; - - // Cf at tracking instant - vector Cfi = Cf00 + m*(Cf - Cf00); - - // Normal vector cross product - vector omega = (n00 ^ n); - - scalar magOmega = mag(omega); - - // magOmega = sin(angle between unit normals) - // Normalise omega vector by magOmega, then multiply by - // angle/dt to give the correct angular velocity vector. - omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue()); - - // Project position onto face and calculate this position - // relative to the face centre. - vector facePos = - p.position() - - ((p.position() - Cfi) & nw)*nw - - Cfi; - - Up += (omega ^ facePos); - } - - // No further action is required if the motion is - // translational only, nw and Up have already been set. } } diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H index 105075ad2f..e980e1cdd5 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -588,8 +588,6 @@ public: ( const parcelType& p, const polyPatch& pp, - const scalar trackFraction, - const tetIndices& tetIs, vector& normal, vector& Up ) const; diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 49bb969210..161b3bb318 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -296,17 +296,15 @@ bool Foam::KinematicParcel::move } else { - // Abandon the track, and move to the end of the sub-step. If the - // the mesh is moving, this will implicitly move the parcel. - if (mesh.moving()) - { - WarningInFunction - << "Tracking was abandoned on a moving mesh. Parcels may " - << "move unphysically as a result." << endl; - } + // At present the only thing that sets active_ to false is a stick + // wall interaction. We want the position of the particle to remain + // the same relative to the face that it is on. The local + // coordinates therefore do not change. We still advance in time and + // perform the relevant interactions with the fixed particle. p.stepFraction() += f; } + const scalar dt = (p.stepFraction() - sfrac)*trackTime; // Avoid problems with extremely small timesteps @@ -325,7 +323,7 @@ bool Foam::KinematicParcel::move if (p.onBoundaryFace() && td.keepParticle) { - if (isA(pbMesh[p.patch(p.face())])) + if (isA(pbMesh[p.patch()])) { td.switchProcessor = true; } diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C index 00393be188..6a15b92a92 100644 --- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C +++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -177,7 +177,7 @@ void Foam::ParticleErosion::postPatch vector Up; // patch-normal direction - this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up); + this->owner().patchData(p, pp, nw, Up); // particle velocity reletive to patch const vector& U = p.U() - Up; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C index 90acd22ca0..a73fdf4a57 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C @@ -237,7 +237,7 @@ bool Foam::LocalInteraction::correct vector nw; vector Up; - this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up); + this->owner().patchData(p, pp, nw, Up); // Calculate motion relative to patch velocity U -= Up; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/MultiInteraction/MultiInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/MultiInteraction/MultiInteraction.C index d2701fbdc8..0ada7622e5 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/MultiInteraction/MultiInteraction.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/MultiInteraction/MultiInteraction.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -171,7 +171,7 @@ bool Foam::MultiInteraction::correct if (p.face() != origFacei) { origFacei = p.face(); - patchi = p.patch(p.face()); + patchi = p.patch(); // Interaction model has moved particle off wall? if (patchi == -1) diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C index f6c1fda460..2ecc54fc38 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ bool Foam::Rebound::correct vector nw; vector Up; - this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up); + this->owner().patchData(p, pp, nw, Up); // Calculate motion relative to patch velocity U -= Up; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C index d7641ac389..a6e43b0829 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C @@ -145,7 +145,7 @@ bool Foam::StandardWallInteraction::correct vector nw; vector Up; - this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up); + this->owner().patchData(p, pp, nw, Up); // Calculate motion relative to patch velocity U -= Up; diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C index d21041a4ec..4b4b8034b9 100644 --- a/src/lagrangian/solidParticle/solidParticle.C +++ b/src/lagrangian/solidParticle/solidParticle.C @@ -84,7 +84,7 @@ bool Foam::solidParticle::move if (onBoundaryFace() && td.keepParticle) { - if (isA(pbMesh[patch(face())])) + if (isA(pbMesh[patch()])) { td.switchProcessor = true; }