diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files index ba19ff80a06e3d3c1aceb2b7ea0e0aee97ebc708..ab0c44e69d1a5716e54a78dea1df3ff9e2a81522 100644 --- a/src/postProcessing/functionObjects/field/Make/files +++ b/src/postProcessing/functionObjects/field/Make/files @@ -33,6 +33,7 @@ wallBoundedStreamLine/wallBoundedStreamLine.C wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C wallBoundedStreamLine/wallBoundedStreamLineParticle.C wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C +wallBoundedStreamLine/wallBoundedParticle.C surfaceInterpolateFields/surfaceInterpolateFields.C surfaceInterpolateFields/surfaceInterpolateFieldsFunctionObject.C diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C index 0cdb3664f9c5e86f7b3ffe23f17b4f999a17afa9..272d38ce59fd070ae044f24b4bbd844f30e78b41 100644 --- a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C +++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C @@ -25,764 +25,9 @@ License #include "wallBoundedStreamLineParticle.H" #include "vectorFieldIOField.H" -#include "polyMesh.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ -// defineParticleTypeNameAndDebug(wallBoundedStreamLineParticle, 0); -} - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -//// Check position is inside tet -//void Foam::wallBoundedStreamLineParticle::checkInside() const -//{ -// const tetIndices ti(currentTetIndices()); -// const tetPointRef tpr(ti.tet(mesh_)); -// if (!tpr.inside(position())) -// { -// FatalErrorIn("wallBoundedStreamLineParticle::checkInside(..)") -// << "Particle:" //<< static_cast<const particle&>(*this) -// << info() -// << "is not inside " << tpr -// << abort(FatalError); -// } -//} -// -// -//void Foam::wallBoundedStreamLineParticle::checkOnEdge() const -//{ -// // Check that edge (as indicated by meshEdgeStart_, diagEdge_) is -// // indeed one that contains the position. -// const edge e = currentEdge(); -// -// linePointRef ln(e.line(mesh_.points())); -// -// pointHit ph(ln.nearestDist(position())); -// -// if (ph.distance() > 1E-6) -// { -// FatalErrorIn -// ( -// "wallBoundedStreamLineParticle::checkOnEdge()" -// ) << "Problem :" -// << " particle:" //<< static_cast<const particle&>(*this) -// << info() -// << "edge:" << e -// << " at:" << ln -// << " distance:" << ph.distance() -// << abort(FatalError); -// } -//} -// -// -//void Foam::wallBoundedStreamLineParticle::checkOnTriangle(const point& p) -//const -//{ -// const triFace tri(currentTetIndices().faceTriIs(mesh_)); -// pointHit ph = tri.nearestPoint(p, mesh_.points()); -// if (ph.distance() > 1e-9) -// { -// FatalErrorIn -// ( -// "wallBoundedStreamLineParticle::checkOnTriangle(const point&)" -// ) << "Problem :" -// << " particle:" //<< static_cast<const particle&>(*this) -// << info() -// << "point:" << p -// << " distance:" << ph.distance() -// << abort(FatalError); -// } -//} - - -// Construct the edge the particle is on (according to meshEdgeStart_, -// diagEdge_) -Foam::edge Foam::wallBoundedStreamLineParticle::currentEdge() const -{ - if ((meshEdgeStart_ != -1) == (diagEdge_ != -1)) - { - FatalErrorIn("wallBoundedStreamLineParticle::currentEdge() const") - << "Particle:" //<< static_cast<const particle&>(*this) - << info() - << "cannot both be on a mesh edge and a face-diagonal edge." - << " meshEdgeStart_:" << meshEdgeStart_ - << " diagEdge_:" << diagEdge_ - << abort(FatalError); - } - - const Foam::face& f = mesh_.faces()[tetFace()]; - - if (meshEdgeStart_ != -1) - { - return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_)); - } - else - { - label faceBasePtI = mesh_.tetBasePtIs()[tetFace()]; - label diagPtI = (faceBasePtI+diagEdge_)%f.size(); - return edge(f[faceBasePtI], f[diagPtI]); - } -} - - -void Foam::wallBoundedStreamLineParticle::crossEdgeConnectedFace -( - const edge& meshEdge -) -{ - //label oldFaceI = tetFace(); - - // Update tetFace, tetPt - particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge); - - // Update face to be same as tracking one - face() = tetFace(); - - - // And adapt meshEdgeStart_. - const Foam::face& f = mesh_.faces()[tetFace()]; - label fp = findIndex(f, meshEdge[0]); - - if (f.nextLabel(fp) == meshEdge[1]) - { - meshEdgeStart_ = fp; - } - else - { - label fpMin1 = f.rcIndex(fp); - - if (f[fpMin1] == meshEdge[1]) - { - meshEdgeStart_ = fpMin1; - } - else - { - FatalErrorIn - ( - "wallBoundedStreamLineParticle::crossEdgeConnectedFace" - "(const edge&)" - ) << "Problem :" - << " particle:" //<< static_cast<const particle&>(*this) - << info() - << "face:" << tetFace() - << " verts:" << f - << " meshEdge:" << meshEdge - << abort(FatalError); - } - } - - diagEdge_ = -1; - - //Pout<< " crossed meshEdge " - // << meshEdge.line(mesh().points()) - // << " from face:" << oldFaceI - // << " to face:" << tetFace() << endl; - - - // Check that still on same mesh edge - - const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_)); - if (eNew != meshEdge) - { - FatalErrorIn - ( - "wallBoundedStreamLineParticle::crossEdgeConnectedFace" - "(const edge&)" - ) << "Problem" << abort(FatalError); - } - - - // Check that edge (as indicated by meshEdgeStart_) is indeed one that - // contains the position. - //checkOnEdge(); -} - - -void Foam::wallBoundedStreamLineParticle::crossDiagonalEdge() -{ - if (diagEdge_ == -1) - { - FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()") - << "Particle:" //<< static_cast<const particle&>(*this) - << info() - << "not on a diagonal edge" << abort(FatalError); - } - if (meshEdgeStart_ != -1) - { - FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()") - << "Particle:" //<< static_cast<const particle&>(*this) - << info() - << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError); - } - - //label oldTetPt = tetPt(); - - const Foam::face& f = mesh_.faces()[tetFace()]; - - // tetPtI starts from 1, goes up to f.size()-2 - - if (tetPt() == diagEdge_) - { - tetPt() = f.rcIndex(tetPt()); - } - else - { - label nextTetPt = f.fcIndex(tetPt()); - if (diagEdge_ == nextTetPt) - { - tetPt() = nextTetPt; - } - else - { - FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()") - << "Particle:" //<< static_cast<const particle&>(*this) - << info() - << "tetPt:" << tetPt() - << " diagEdge:" << diagEdge_ << abort(FatalError); - } - } - - meshEdgeStart_ = -1; - - //Pout<< " crossed diagonalEdge " - // << currentEdge().line(mesh().points()) - // << " from tetPt:" << oldTetPt - // << " to tetPt:" << tetPt() << endl; -} - - -//- Track through a single triangle. -// Gets passed tet+triangle the particle is in. Updates position() but nothing -// else. Returns the triangle edge the particle is now on. -Foam::scalar Foam::wallBoundedStreamLineParticle::trackFaceTri -( - const vector& endPosition, - label& minEdgeI -) -{ - // Track p from position to endPosition - const triFace tri(currentTetIndices().faceTriIs(mesh_)); - vector n = tri.normal(mesh_.points()); - //if (mag(n) < sqr(SMALL)) - //{ - // FatalErrorIn("wallBoundedStreamLineParticle::trackFaceTri(..)") - // << "Small triangle." //<< static_cast<const particle&>(*this) - // << info() - // << "n:" << n - // << abort(FatalError); - //} - n /= mag(n)+VSMALL; - - // Check which edge intersects the trajectory. - // Project trajectory onto triangle - minEdgeI = -1; - scalar minS = 1; // end position - - //const point oldPosition(position()); - - - edge currentE(-1, -1); - if (meshEdgeStart_ != -1 || diagEdge_ != -1) - { - currentE = currentEdge(); - } - - // Determine path along line position+s*d to see where intersections - // are. - - forAll(tri, i) - { - label j = tri.fcIndex(i); - - const point& pt0 = mesh_.points()[tri[i]]; - const point& pt1 = mesh_.points()[tri[j]]; - - if (edge(tri[i], tri[j]) == currentE) - { - // Do not check particle is on - continue; - } - - // Outwards pointing normal - vector edgeNormal = (pt1-pt0)^n; - - //if (mag(edgeNormal) < SMALL) - //{ - // FatalErrorIn("wallBoundedStreamLineParticle::trackFaceTri(..)") - // << "Edge not perpendicular to triangle." - // //<< static_cast<const particle&>(*this) - // << info() - // << "triangle n:" << n - // << " edgeNormal:" << edgeNormal - // << " on tri:" << tri - // << " at:" << pt0 - // << " at:" << pt1 - // << abort(FatalError); - //} - - - edgeNormal /= mag(edgeNormal)+VSMALL; - - // Determine whether position and end point on either side of edge. - scalar sEnd = (endPosition-pt0)&edgeNormal; - if (sEnd >= 0) - { - // endPos is outside triangle. position() should always be - // inside. - scalar sStart = (position()-pt0)&edgeNormal; - if (mag(sEnd-sStart) > VSMALL) - { - scalar s = sStart/(sStart-sEnd); - - if (s >= 0 && s < minS) - { - minS = s; - minEdgeI = i; - } - } - } - } - - if (minEdgeI != -1) - { - position() += minS*(endPosition-position()); - } - else - { - // Did not hit any edge so tracked to the end position. Set position - // without any calculation to avoid truncation errors. - position() = endPosition; - minS = 1.0; - } - - // Project position() back onto plane of triangle - const point& triPt = mesh_.points()[tri[0]]; - position() -= ((position()-triPt)&n)*n; - - - //Pout<< " tracked from:" << oldPosition << " to:" << position() - // << " projectedEnd:" << endPosition - // << " at s:" << minS << endl; - //if (minEdgeI != -1) - //{ - // Pout<< " on edge:" << minEdgeI - // << " on edge:" - // << mesh_.points()[tri[minEdgeI]] - // << mesh_.points()[tri[tri.fcIndex(minEdgeI)]] - // << endl; - //} - - return minS; -} - - -// See if the current triangle has got a point on the -// correct side of the edge. -bool Foam::wallBoundedStreamLineParticle::isTriAlongTrack -( - const point& endPosition -) const -{ - const triFace triVerts(currentTetIndices().faceTriIs(mesh_)); - const edge currentE = currentEdge(); - - //if (debug) - { - if - ( - currentE[0] == currentE[1] - || findIndex(triVerts, currentE[0]) == -1 - || findIndex(triVerts, currentE[1]) == -1 - ) - { - FatalErrorIn - ( - "wallBoundedStreamLineParticle::isTriAlongTrack" - "(const point&)" - ) << "Edge " << currentE << " not on triangle " << triVerts - << info() - << abort(FatalError); - } - } - - - const vector dir = endPosition-position(); - - // Get normal of currentE - vector n = triVerts.normal(mesh_.points()); - n /= mag(n); - - forAll(triVerts, i) - { - label j = triVerts.fcIndex(i); - const point& pt0 = mesh_.points()[triVerts[i]]; - const point& pt1 = mesh_.points()[triVerts[j]]; - - if (edge(triVerts[i], triVerts[j]) == currentE) - { - vector edgeNormal = (pt1-pt0)^n; - return (dir&edgeNormal) < 0; - } - } - - FatalErrorIn - ( - "wallBoundedStreamLineParticle::isTriAlongTrack" - "(const point&)" - ) << "Problem" << abort(FatalError); - - return false; -} - - -void Foam::wallBoundedStreamLineParticle::patchInteraction -( - trackingData& td, - const scalar trackFraction -) -{ - typedef trackingData::cloudType cloudType; - typedef cloudType::particleType particleType; - - particleType& p = static_cast<particleType&>(*this); - p.hitFace(td); - - if (!internalFace(faceI_)) - { - label origFaceI = faceI_; - label patchI = patch(faceI_); - - // No action taken for tetPtI_ for tetFaceI_ here, handled by - // patch interaction call or later during processor transfer. - - - // Dummy tet indices. What to do here? - tetIndices faceHitTetIs; - - if - ( - !p.hitPatch - ( - mesh_.boundaryMesh()[patchI], - td, - patchI, - trackFraction, - faceHitTetIs - ) - ) - { - // Did patch interaction model switch patches? - // Note: recalculate meshEdgeStart_, diagEdge_! - if (faceI_ != origFaceI) - { - patchI = patch(faceI_); - } - - const polyPatch& patch = mesh_.boundaryMesh()[patchI]; - - if (isA<wedgePolyPatch>(patch)) - { - p.hitWedgePatch - ( - static_cast<const wedgePolyPatch&>(patch), td - ); - } - else if (isA<symmetryPolyPatch>(patch)) - { - p.hitSymmetryPatch - ( - static_cast<const symmetryPolyPatch&>(patch), td - ); - } - else if (isA<cyclicPolyPatch>(patch)) - { - p.hitCyclicPatch - ( - static_cast<const cyclicPolyPatch&>(patch), td - ); - } - else if (isA<processorPolyPatch>(patch)) - { - p.hitProcessorPatch - ( - static_cast<const processorPolyPatch&>(patch), td - ); - } - else if (isA<wallPolyPatch>(patch)) - { - p.hitWallPatch - ( - static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs - ); - } - else - { - p.hitPatch(patch, td); - } - } - } -} - -//- Track particle to a given position and returns 1.0 if the -// trajectory is completed without hitting a face otherwise -// stops at the face and returns the fraction of the trajectory -// completed. -// on entry 'stepFraction()' should be set to the fraction of the -// time-step at which the tracking starts. -Foam::scalar Foam::wallBoundedStreamLineParticle::trackToEdge -( - trackingData& td, - const vector& endPosition -) -{ - // Are we on a track face? If not we do a topological walk. - - // Particle: - // - cell_ always set - // - tetFace_, tetPt_ always set (these identify tet particle is in) - // - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on) - - //checkInside(); - //checkOnTriangle(position()); - //if (meshEdgeStart_ != -1 || diagEdge_ != -1) - //{ - // checkOnEdge(); - //} - - scalar trackFraction = 0.0; - - if (!td.isWallPatch_[tetFace()]) - { - // Don't track across face. Just walk in cell. Particle is on - // mesh edge (as indicated by meshEdgeStart_). - const edge meshEdge(currentEdge()); - - // If internal face check whether to go to neighbour cell or just - // check to the other internal tet on the edge. - if (mesh_.isInternalFace(tetFace())) - { - label nbrCellI = - ( - cellI_ == mesh_.faceOwner()[faceI_] - ? mesh_.faceNeighbour()[faceI_] - : mesh_.faceOwner()[faceI_] - ); - // Check angle to nbrCell tet. Is it in the direction of the - // endposition? I.e. since volume of nbr tet is positive the - // tracking direction should be into the tet. - tetIndices nbrTi(nbrCellI, tetFaceI_, tetPtI_, mesh_); - if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0) - { - // Change into nbrCell. No need to change tetFace, tetPt. - //Pout<< " crossed from cell:" << cellI_ - // << " into " << nbrCellI << endl; - cellI_ = nbrCellI; - patchInteraction(td, trackFraction); - } - else - { - // Walk to other face on edge. Changes tetFace, tetPt but not - // cell. - crossEdgeConnectedFace(meshEdge); - patchInteraction(td, trackFraction); - } - } - else - { - // Walk to other face on edge. This might give loop since - // particle should have been removed? - crossEdgeConnectedFace(meshEdge); - patchInteraction(td, trackFraction); - } - } - else - { - // We're inside a tet on the wall. Check if the current tet is - // the one to cross. If not we cross into the neighbouring triangle. - - if (mesh_.isInternalFace(tetFace())) - { - FatalErrorIn - ( - "wallBoundedStreamLineParticle::trackToEdge" - "(trackingData&, const vector&)" - ) << "Can only track on boundary faces." - << " Face:" << tetFace() - << " at:" << mesh_.faceCentres()[tetFace()] - << abort(FatalError); - } - - point projectedEndPosition = endPosition; - // Remove normal component - { - const triFace tri(currentTetIndices().faceTriIs(mesh_)); - vector n = tri.normal(mesh_.points()); - n /= mag(n); - const point& basePt = mesh_.points()[tri[0]]; - projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n; - } - - - bool doTrack = false; - if (meshEdgeStart_ == -1 && diagEdge_ == -1) - { - // We're starting and not yet on an edge. - doTrack = true; - } - else - { - // See if the current triangle has got a point on the - // correct side of the edge. - doTrack = isTriAlongTrack(projectedEndPosition); - } - - - if (doTrack) - { - // Track across triangle. Return triangle edge crossed. - label triEdgeI = -1; - trackFraction = trackFaceTri(projectedEndPosition, triEdgeI); - - if (triEdgeI == -1) - { - // Reached endpoint - //checkInside(); - diagEdge_ = -1; - meshEdgeStart_ = -1; - return trackFraction; - } - - const tetIndices ti(currentTetIndices()); - - // Triangle (faceTriIs) gets constructed from - // f[faceBasePtI_], - // f[facePtAI_], - // f[facePtBI_] - // - // So edge indices are: - // 0 : edge between faceBasePtI_ and facePtAI_ - // 1 : edge between facePtAI_ and facePtBI_ (is always a real edge) - // 2 : edge between facePtBI_ and faceBasePtI_ - - const Foam::face& f = mesh_.faces()[ti.face()]; - const label fp0 = ti.faceBasePt(); - - if (triEdgeI == 0) - { - if (ti.facePtA() == f.fcIndex(fp0)) - { - //Pout<< "Real edge." << endl; - diagEdge_ = -1; - meshEdgeStart_ = fp0; - //checkOnEdge(); - crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); - } - else if (ti.facePtA() == f.rcIndex(fp0)) - { - //Note: should not happen since boundary face so owner - //Pout<< "Real edge." << endl; - FatalErrorIn("shold not happend") << info() - << abort(FatalError); - - diagEdge_ = -1; - meshEdgeStart_ = f.rcIndex(fp0); - //checkOnEdge(); - crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); - } - else - { - // Get index of triangle on other side of edge. - diagEdge_ = ti.facePtA()-fp0; - if (diagEdge_ < 0) - { - diagEdge_ += f.size(); - } - meshEdgeStart_ = -1; - //checkOnEdge(); - crossDiagonalEdge(); - } - } - else if (triEdgeI == 1) - { - //Pout<< "Real edge." << endl; - diagEdge_ = -1; - meshEdgeStart_ = ti.facePtA(); - //checkOnEdge(); - crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); - } - else // if (triEdgeI == 2) - { - if (ti.facePtB() == f.rcIndex(fp0)) - { - //Pout<< "Real edge." << endl; - diagEdge_ = -1; - meshEdgeStart_ = ti.facePtB(); - //checkOnEdge(); - crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); - } - else if (ti.facePtB() == f.fcIndex(fp0)) - { - //Note: should not happen since boundary face so owner - //Pout<< "Real edge." << endl; - FatalErrorIn("shold not happend") << info() - << abort(FatalError); - - diagEdge_ = -1; - meshEdgeStart_ = fp0; - //checkOnEdge(); - crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); - } - else - { - //Pout<< "Triangle edge." << endl; - // Get index of triangle on other side of edge. - diagEdge_ = ti.facePtB()-fp0; - if (diagEdge_ < 0) - { - diagEdge_ += f.size(); - } - meshEdgeStart_ = -1; - //checkOnEdge(); - crossDiagonalEdge(); - } - } - } - else - { - // Current tet is not the right one. Check the neighbour tet. - - if (meshEdgeStart_ != -1) - { - // Particle is on mesh edge so change into other face on cell - crossEdgeConnectedFace(currentEdge()); - //checkOnEdge(); - patchInteraction(td, trackFraction); - } - else - { - // Particle is on diagonal edge so change into the other - // triangle. - crossDiagonalEdge(); - //checkOnEdge(); - } - } - } - - //checkInside(); - - return trackFraction; -} - - Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields ( const trackingData& td, @@ -897,22 +142,18 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle const label lifeTime ) : - particle(mesh, position, cellI, tetFaceI, tetPtI), - meshEdgeStart_(meshEdgeStart), - diagEdge_(diagEdge), + wallBoundedParticle + ( + mesh, + position, + cellI, + tetFaceI, + tetPtI, + meshEdgeStart, + diagEdge + ), lifeTime_(lifeTime) -{ - //checkInside(); - - //if (meshEdgeStart_ != -1 || diagEdge_ != -1) - //{ - // checkOnEdge(); - //} - - // Unfortunately have no access to trackdata so cannot check if particle - // is on a wallPatch or has an mesh edge set (either of which is - // a requirement). -} +{} Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle @@ -922,14 +163,14 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle bool readFields ) : - particle(mesh, is, readFields) + wallBoundedParticle(mesh, is, readFields) { if (readFields) { List<scalarList> sampledScalars; List<vectorList> sampledVectors; - is >> meshEdgeStart_ >> diagEdge_ >> lifeTime_ + is >> lifeTime_ >> sampledPositions_ >> sampledScalars >> sampledVectors; sampledScalars_.setSize(sampledScalars.size()); @@ -958,9 +199,7 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle const wallBoundedStreamLineParticle& p ) : - particle(p), - meshEdgeStart_(p.meshEdgeStart_), - diagEdge_(p.diagEdge_), + wallBoundedParticle(p), lifeTime_(p.lifeTime_), sampledPositions_(p.sampledPositions_), sampledScalars_(p.sampledScalars_), @@ -1102,101 +341,6 @@ bool Foam::wallBoundedStreamLineParticle::move } -bool Foam::wallBoundedStreamLineParticle::hitPatch -( - const polyPatch&, - trackingData& td, - const label patchI, - const scalar trackFraction, - const tetIndices& tetIs -) -{ - // Disable generic patch interaction - return false; -} - - -void Foam::wallBoundedStreamLineParticle::hitWedgePatch -( - const wedgePolyPatch& pp, - trackingData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -void Foam::wallBoundedStreamLineParticle::hitSymmetryPatch -( - const symmetryPolyPatch& pp, - trackingData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -void Foam::wallBoundedStreamLineParticle::hitCyclicPatch -( - const cyclicPolyPatch& pp, - trackingData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -void Foam::wallBoundedStreamLineParticle::hitProcessorPatch -( - const processorPolyPatch& pp, - trackingData& td -) -{ - // Switch particle - td.switchProcessor = true; - - // Adapt edgeStart_ for other side. - // E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so - // on the other side between 2 and 3 so edgeStart_ should be - // f.size()-edgeStart_-1. - - const Foam::face& f = mesh_.faces()[face()]; - - if (meshEdgeStart_ != -1) - { - meshEdgeStart_ = f.size()-meshEdgeStart_-1; - } - else - { - // diagEdge_ is relative to faceBasePt - diagEdge_ = f.size()-diagEdge_; - } -} - - -void Foam::wallBoundedStreamLineParticle::hitWallPatch -( - const wallPolyPatch& wpp, - trackingData& td, - const tetIndices& -) -{} - - -void Foam::wallBoundedStreamLineParticle::hitPatch -( - const polyPatch& wpp, - trackingData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - void Foam::wallBoundedStreamLineParticle::readFields ( Cloud<wallBoundedStreamLineParticle>& c @@ -1207,18 +351,7 @@ void Foam::wallBoundedStreamLineParticle::readFields return; } - particle::readFields(c); - - IOField<label> meshEdgeStart - ( - c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ) - ); - - IOField<label> diagEdge - ( - c.fieldIOobject("diagEdge_", IOobject::MUST_READ) - ); - c.checkFieldIOobject(c, diagEdge); + wallBoundedParticle::readFields(c); IOField<label> lifeTime ( @@ -1235,8 +368,6 @@ void Foam::wallBoundedStreamLineParticle::readFields label i = 0; forAllIter(Cloud<wallBoundedStreamLineParticle>, c, iter) { - iter().meshEdgeStart_ = meshEdgeStart[i]; - iter().diagEdge_ = diagEdge[i]; iter().lifeTime_ = lifeTime[i]; iter().sampledPositions_.transfer(sampledPositions[i]); i++; @@ -1249,20 +380,10 @@ void Foam::wallBoundedStreamLineParticle::writeFields const Cloud<wallBoundedStreamLineParticle>& c ) { - particle::writeFields(c); + wallBoundedParticle::writeFields(c); label np = c.size(); - IOField<label> meshEdgeStart - ( - c.fieldIOobject("meshEdgeStart", IOobject::NO_READ), - np - ); - IOField<label> diagEdge - ( - c.fieldIOobject("diagEdge", IOobject::NO_READ), - np - ); IOField<label> lifeTime ( c.fieldIOobject("lifeTime", IOobject::NO_READ), @@ -1277,15 +398,11 @@ void Foam::wallBoundedStreamLineParticle::writeFields label i = 0; forAllConstIter(Cloud<wallBoundedStreamLineParticle>, c, iter) { - meshEdgeStart[i] = iter().meshEdgeStart_; - diagEdge[i] = iter().diagEdge_; lifeTime[i] = iter().lifeTime_; sampledPositions[i] = iter().sampledPositions_; i++; } - meshEdgeStart.write(); - diagEdge.write(); lifeTime.write(); sampledPositions.write(); } @@ -1299,9 +416,7 @@ Foam::Ostream& Foam::operator<< const wallBoundedStreamLineParticle& p ) { - os << static_cast<const particle&>(p) - << token::SPACE << p.meshEdgeStart_ - << token::SPACE << p.diagEdge_ + os << static_cast<const wallBoundedParticle&>(p) << token::SPACE << p.lifeTime_ << token::SPACE << p.sampledPositions_ << token::SPACE << p.sampledScalars_ @@ -1317,38 +432,4 @@ Foam::Ostream& Foam::operator<< } -Foam::Ostream& Foam::operator<< -( - Ostream& os, - const InfoProxy<wallBoundedStreamLineParticle>& ip -) -{ - const wallBoundedStreamLineParticle& p = ip.t_; - - tetPointRef tpr(p.currentTet()); - - os << " " << static_cast<const particle&>(p) << nl - << " tet:" << nl; - os << " "; - meshTools::writeOBJ(os, tpr.a()); - os << " "; - meshTools::writeOBJ(os, tpr.b()); - os << " "; - meshTools::writeOBJ(os, tpr.c()); - os << " "; - meshTools::writeOBJ(os, tpr.d()); - os << " l 1 2" << nl - << " l 1 3" << nl - << " l 1 4" << nl - << " l 2 3" << nl - << " l 2 4" << nl - << " l 3 4" << nl; - os << " "; - meshTools::writeOBJ(os, p.position()); - - return os; -} - - - // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H index a70311c4cfc592250a4ea02a1fe1340485f9a005..54e6a4760adc177652a0cdea4c1e4230aa77d201 100644 --- a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H +++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H @@ -36,7 +36,7 @@ SourceFiles #ifndef wallBoundedStreamLineParticle_H #define wallBoundedStreamLineParticle_H -#include "particle.H" +#include "wallBoundedParticle.H" #include "autoPtr.H" #include "interpolation.H" #include "vectorList.H" @@ -50,20 +50,23 @@ namespace Foam class wallBoundedStreamLineParticleCloud; /*---------------------------------------------------------------------------*\ - Class wallBoundedStreamLineParticle Declaration + Class wallBoundedStreamLineParticle Declaration \*---------------------------------------------------------------------------*/ class wallBoundedStreamLineParticle : - public particle + public wallBoundedParticle { public: - //- Class used to pass tracking data to the trackToFace function + //- Class used to pass tracking data to the trackToEdge function class trackingData : - public particle::TrackingData<Cloud<wallBoundedStreamLineParticle> > + public wallBoundedParticle::TrackingData + < + Cloud<wallBoundedStreamLineParticle> + > { public: @@ -75,8 +78,6 @@ public: const bool trackForward_; const scalar trackLength_; - const PackedBoolList& isWallPatch_; - DynamicList<vectorList>& allPositions_; List<DynamicList<scalarList> >& allScalars_; List<DynamicList<vectorList> >& allVectors_; @@ -99,16 +100,19 @@ public: List<DynamicList<vectorList> >& allVectors ) : - particle::TrackingData<Cloud<wallBoundedStreamLineParticle> > + wallBoundedParticle::TrackingData + < + Cloud<wallBoundedStreamLineParticle> + > ( - cloud + cloud, + isWallPatch ), vsInterp_(vsInterp), vvInterp_(vvInterp), UIndex_(UIndex), trackForward_(trackForward), trackLength_(trackLength), - isWallPatch_(isWallPatch), allPositions_(allPositions), allScalars_(allScalars), @@ -121,21 +125,6 @@ private: // Private data - //- Particle is on mesh edge: - // const face& f = mesh.faces()[tetFace()] - // const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_)); - // Note that this real edge - // is also one of the edges of the face-triangle (from - // tetFace()+tetPt()). - label meshEdgeStart_; - - //- Particle is on diagonal edge: - // const face& f = mesh.faces()[tetFace()] - // label faceBasePtI = mesh.tetBasePtIs()[faceI]; - // label diagPtI = (faceBasePtI+diagEdge_)%f.size(); - // const edge e(f[faceBasePtI], f[diagPtI]); - label diagEdge_; - //- Lifetime of particle. Particle dies when reaches 0. label lifeTime_; @@ -151,52 +140,17 @@ private: // Private Member Functions - //- Construct current edge - edge currentEdge() const; - - //- Check if inside current tet - //void checkInside() const; - - //- Check if on current edge - //void checkOnEdge() const; - - //- Check if point on triangle - //void checkOnTriangle(const point&) const; - - //- Cross mesh edge into different face on same cell - void crossEdgeConnectedFace(const edge& meshEdge); - - //- Cross diagonal edge into different triangle on same face,cell - void crossDiagonalEdge(); - - //- Track through single triangle - scalar trackFaceTri(const vector& endPosition, label& minEdgeI); - - //- Do all patch interaction - void patchInteraction(trackingData& td, const scalar trackFraction); - - //- Is current triangle in the track direction - bool isTriAlongTrack(const point& endPosition) const; - - //- Equivalent of trackToFace - scalar trackToEdge - ( - trackingData& td, - const vector& endPosition - ); - - //- Interpolate all quantities; return interpolated velocity. vector interpolateFields ( - const trackingData&, - const point&, + const trackingData& td, + const point& position, const label cellI, const label faceI ); - //- Interpolate all quantities, return updated velocity. vector sample(trackingData& td); + public: // Constructors @@ -265,73 +219,6 @@ public: bool move(trackingData&, const scalar trackTime); - //- Overridable function to handle the particle hitting a patch - // Executed before other patch-hitting functions - bool hitPatch - ( - const polyPatch&, - trackingData& td, - const label patchI, - const scalar trackFraction, - const tetIndices& tetIs - ); - - //- Overridable function to handle the particle hitting a wedge - void hitWedgePatch - ( - const wedgePolyPatch&, - trackingData& td - ); - - //- Overridable function to handle the particle hitting a - // symmetryPlane - void hitSymmetryPatch - ( - const symmetryPolyPatch&, - trackingData& td - ); - - //- Overridable function to handle the particle hitting a cyclic - void hitCyclicPatch - ( - const cyclicPolyPatch&, - trackingData& td - ); - - //- Overridable function to handle the particle hitting a - //- processorPatch - void hitProcessorPatch - ( - const processorPolyPatch&, - trackingData& td - ); - - //- Overridable function to handle the particle hitting a wallPatch - void hitWallPatch - ( - const wallPolyPatch&, - trackingData& td, - const tetIndices& - ); - - //- Overridable function to handle the particle hitting a polyPatch - void hitPatch - ( - const polyPatch&, - trackingData& td - ); - - - // Info - - //- Return info proxy. - // Used to print particle information to a stream - InfoProxy<wallBoundedStreamLineParticle> info() const - { - return *this; - } - - // I-O //- Read @@ -351,12 +238,6 @@ public: Ostream&, const wallBoundedStreamLineParticle& ); - - friend Ostream& operator<< - ( - Ostream&, - const InfoProxy<wallBoundedStreamLineParticle>& - ); };