diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H index 8ab5d1e2e8..0ddf19294c 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H @@ -142,18 +142,30 @@ public: inline label& tetPt(); //- Return the indices corresponding to the tri on the face for - // this tet. The normal of the tri points out of the cell. - inline triFace faceTriIs(const polyMesh& mesh) const; + // this tet. The normal of the tri points out of the cell + inline triFace faceTriIs + ( + const polyMesh& mesh, + const bool warn = true + ) const; + + //- Return the local indices corresponding to the tri on the face + // for this tet. The normal of the tri points out of the cell + inline triFace triIs + ( + const polyMesh& mesh, + const bool warn = true + ) const; //- Return the geometry corresponding to this tet inline tetPointRef tet(const polyMesh& mesh) const; //- Return the geometry corresponding to the tri on the face for - // this tet. The normal of the tri points out of the cell. + // this tet. The normal of the tri points out of the cell inline triPointRef faceTri(const polyMesh& mesh) const; //- Return the geometry corresponding to the tri on the face for - // this tet using the old positions. + // this tet using the old positions inline triPointRef oldFaceTri(const polyMesh& mesh) const; diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H index 70edef09b1..d09eeeefea 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H @@ -61,7 +61,11 @@ inline Foam::label& Foam::tetIndices::tetPt() } -inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const +inline Foam::triFace Foam::tetIndices::faceTriIs +( + const polyMesh& mesh, + const bool warn +) const { const Foam::face& f = mesh.faces()[face()]; @@ -70,22 +74,26 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const if (faceBasePtI < 0) { faceBasePtI = 0; - if (nWarnings < maxNWarnings) + + if (warn) { - WarningInFunction - << "No base point for face " << face() << ", " << f - << ", produces a valid tet decomposition." << endl; - ++ nWarnings; - } - if (nWarnings == maxNWarnings) - { - Warning - << "Suppressing any further warnings." << endl; - ++ nWarnings; + if (nWarnings < maxNWarnings) + { + WarningInFunction + << "No base point for face " << face() << ", " << f + << ", produces a valid tet decomposition." << endl; + ++nWarnings; + } + if (nWarnings == maxNWarnings) + { + Warning + << "Suppressing any further warnings." << endl; + ++nWarnings; + } } } - label facePtI = (tetPt() + faceBasePtI) % f.size(); + label facePtI = (tetPti_ + faceBasePtI) % f.size(); label faceOtherPtI = f.fcIndex(facePtI); if (mesh.faceOwner()[face()] != cell()) @@ -97,6 +105,50 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const } +inline Foam::triFace Foam::tetIndices::triIs +( + const polyMesh& mesh, + const bool warn +) const +{ + const Foam::face& f = mesh.faces()[face()]; + + label faceBasePtI = mesh.tetBasePtIs()[face()]; + + if (faceBasePtI < 0) + { + faceBasePtI = 0; + + if (warn) + { + if (nWarnings < maxNWarnings) + { + WarningInFunction + << "No base point for face " << face() << ", " << f + << ", produces a valid tet decomposition." << endl; + ++nWarnings; + } + if (nWarnings == maxNWarnings) + { + Warning + << "Suppressing any further warnings." << endl; + ++nWarnings; + } + } + } + + label facePtI = (tetPti_ + faceBasePtI) % f.size(); + label faceOtherPtI = f.fcIndex(facePtI); + + if (mesh.faceOwner()[face()] != cell()) + { + Swap(facePtI, faceOtherPtI); + } + + return triFace(faceBasePtI, facePtI, faceOtherPtI); +} + + inline Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const { const pointField& meshPoints = mesh.points(); diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files index 3362b2dbca..b44dc3819f 100644 --- a/src/functionObjects/field/Make/files +++ b/src/functionObjects/field/Make/files @@ -25,12 +25,10 @@ streamLine/streamLineBase.C streamLine/streamLineParticle.C streamLine/streamLineParticleCloud.C -/* wallBoundedStreamLine/wallBoundedStreamLine.C wallBoundedStreamLine/wallBoundedStreamLineParticle.C wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C wallBoundedStreamLine/wallBoundedParticle.C -*/ surfaceInterpolate/surfaceInterpolate.C diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C index 681956f30f..03c1aea82b 100644 --- a/src/functionObjects/field/nearWallFields/nearWallFields.C +++ b/src/functionObjects/field/nearWallFields/nearWallFields.C @@ -111,11 +111,12 @@ void Foam::functionObjects::nearWallFields::calcAddressing() tetPti = (startInfo.index()+1) % f.size(); start = startInfo.hitPoint(); + //// Uncomment below to shift slightly in: - //tetIndices tet(celli, meshFacei, tetPti, mesh_); - //start = - // (1.0-1e-6)*startInfo.hitPoint() - // + 1e-6*tet.tet(mesh_).centre(); + tetIndices tet(celli, meshFacei, tetPti); + start = + (1.0 - 1e-6)*startInfo.hitPoint() + + 1e-6*tet.tet(mesh_).centre(); } else { diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C index 6996e1500c..8a21bba822 100644 --- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C +++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C @@ -35,57 +35,6 @@ const std::size_t Foam::wallBoundedParticle::sizeofFields_ // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // -Foam::tetIndices Foam::wallBoundedParticle::currentTetIndices() const -{ - // Replacement for particle::currentTetIndices that avoids error - // upon invalid tetBasePtIs - - const faceList& pFaces = mesh_.faces(); - const labelList& pOwner = mesh_.faceOwner(); - - const Foam::face& f = pFaces[tetFacei_]; - - bool own = (pOwner[tetFacei_] == celli_); - - label faceBasePtI = mesh_.tetBasePtIs()[tetFacei_]; - if (faceBasePtI == -1) - { - //WarningInFunction - // << "No base point for face " << tetFacei_ << ", " << f - // << ", produces a decomposition that has a minimum " - // << "volume greater than tolerance." - // << endl; - faceBasePtI = 0; - } - - label facePtI = (tetPti_ + faceBasePtI) % f.size(); - label otherFacePtI = f.fcIndex(facePtI); - - label facePtAI; - label facePtBI; - if (own) - { - facePtAI = facePtI; - facePtBI = otherFacePtI; - } - else - { - facePtAI = otherFacePtI; - facePtBI = facePtI; - } - - return tetIndices - ( - celli_, - tetFacei_, - faceBasePtI, - facePtAI, - facePtBI, - tetPti_ - ); -} - - Foam::edge Foam::wallBoundedParticle::currentEdge() const { if ((meshEdgeStart_ != -1) == (diagEdge_ != -1)) @@ -99,7 +48,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const << abort(FatalError); } - const Foam::face& f = mesh_.faces()[tetFace()]; + const Foam::face& f = mesh().faces()[tetFace()]; if (meshEdgeStart_ != -1) { @@ -107,7 +56,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const } else { - label faceBasePti = mesh_.tetBasePtIs()[tetFace()]; + label faceBasePti = mesh().tetBasePtIs()[tetFace()]; if (faceBasePti == -1) { //FatalErrorInFunction @@ -120,7 +69,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const faceBasePti = 0; } - label diagPti = (faceBasePti+diagEdge_)%f.size(); + label diagPti = (faceBasePti + diagEdge_)%f.size(); return edge(f[faceBasePti], f[diagPti]); } @@ -135,26 +84,24 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace const edge& e ) { - const faceList& pFaces = mesh_.faces(); - const cellList& pCells = mesh_.cells(); + const faceList& pFaces = mesh().faces(); + const cellList& pCells = mesh().cells(); const Foam::face& f = pFaces[tetFacei]; const Foam::cell& thisCell = pCells[celli]; - forAll(thisCell, cFI) + for (const label facei : thisCell) { // Loop over all other faces of this cell and // find the one that shares this edge - label fI = thisCell[cFI]; - - if (tetFacei == fI) + if (tetFacei == facei) { continue; } - const Foam::face& otherFace = pFaces[fI]; + const Foam::face& otherFace = pFaces[facei]; label edDir = otherFace.edgeDirection(e); @@ -162,7 +109,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace { continue; } - else if (f == pFaces[fI]) + else if (f == pFaces[facei]) { // This is a necessary condition if using duplicate baffles // (so coincident faces). We need to make sure we don't cross into @@ -174,7 +121,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace else { //Found edge on other face - tetFacei = fI; + tetFacei = facei; label eIndex = -1; @@ -192,7 +139,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace eIndex = findIndex(otherFace, e.end()); } - label tetBasePtI = mesh_.tetBasePtIs()[fI]; + label tetBasePtI = mesh().tetBasePtIs()[facei]; if (tetBasePtI == -1) { @@ -246,7 +193,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace(const edge& meshEdge) face() = tetFace(); // And adapt meshEdgeStart_. - const Foam::face& f = mesh_.faces()[tetFace()]; + const Foam::face& f = mesh().faces()[tetFace()]; label fp = findIndex(f, meshEdge[0]); if (f.nextLabel(fp) == meshEdge[1]) @@ -303,7 +250,7 @@ void Foam::wallBoundedParticle::crossDiagonalEdge() << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError); } - const Foam::face& f = mesh_.faces()[tetFace()]; + const Foam::face& f = mesh().faces()[tetFace()]; // tetPti starts from 1, goes up to f.size()-2 @@ -340,7 +287,7 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri ) { // Track p from position to endPosition - const triFace tri(currentTetIndices().faceTriIs(mesh_)); + const triFace tri(currentTetIndices().faceTriIs(mesh(), false)); // Check which edge intersects the trajectory. // Project trajectory onto triangle @@ -358,8 +305,8 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri { label j = tri.fcIndex(i); - const point& pt0 = mesh_.points()[tri[i]]; - const point& pt1 = mesh_.points()[tri[j]]; + const point& pt0 = mesh().points()[tri[i]]; + const point& pt1 = mesh().points()[tri[j]]; if (edge(tri[i], tri[j]) == currentE) { @@ -368,20 +315,20 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri } // Outwards pointing normal - vector edgeNormal = (pt1-pt0)^n; + vector edgeNormal = (pt1 - pt0)^n; - edgeNormal /= mag(edgeNormal)+VSMALL; + edgeNormal /= mag(edgeNormal) + VSMALL; // Determine whether position and end point on either side of edge. - scalar sEnd = (endPosition-pt0)&edgeNormal; + scalar sEnd = (endPosition - pt0)&edgeNormal; if (sEnd >= 0) { - // endPos is outside triangle. position() should always be + // endPos is outside triangle. localPosition_ should always be // inside. - scalar sStart = (position()-pt0)&edgeNormal; - if (mag(sEnd-sStart) > VSMALL) + scalar sStart = (localPosition_ - pt0)&edgeNormal; + if (mag(sEnd - sStart) > VSMALL) { - scalar s = sStart/(sStart-sEnd); + scalar s = sStart/(sStart - sEnd); if (s >= 0 && s < minS) { @@ -394,19 +341,19 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri if (minEdgei != -1) { - position() += minS*(endPosition-position()); + localPosition_ += minS*(endPosition - localPosition_); } else { // Did not hit any edge so tracked to the end position. Set position // without any calculation to avoid truncation errors. - position() = endPosition; + localPosition_ = endPosition; minS = 1.0; } - // Project position() back onto plane of triangle - const point& triPt = mesh_.points()[tri[0]]; - position() -= ((position()-triPt)&n)*n; + // Project localPosition_ back onto plane of triangle + const point& triPt = mesh().points()[tri[0]]; + localPosition_ -= ((localPosition_ - triPt)&n)*n; return minS; } @@ -418,7 +365,7 @@ bool Foam::wallBoundedParticle::isTriAlongTrack const point& endPosition ) const { - const triFace triVerts(currentTetIndices().faceTriIs(mesh_)); + const triFace triVerts(currentTetIndices().faceTriIs(mesh(), false)); const edge currentE = currentEdge(); if @@ -435,13 +382,13 @@ bool Foam::wallBoundedParticle::isTriAlongTrack } - const vector dir = endPosition-position(); + const vector dir = endPosition - localPosition_; forAll(triVerts, i) { label j = triVerts.fcIndex(i); - const point& pt0 = mesh_.points()[triVerts[i]]; - const point& pt1 = mesh_.points()[triVerts[j]]; + const point& pt0 = mesh().points()[triVerts[i]]; + const point& pt1 = mesh().points()[triVerts[j]]; if (edge(triVerts[i], triVerts[j]) == currentE) { @@ -462,7 +409,7 @@ bool Foam::wallBoundedParticle::isTriAlongTrack Foam::wallBoundedParticle::wallBoundedParticle ( const polyMesh& mesh, - const vector& position, + const point& position, const label celli, const label tetFacei, const label tetPti, @@ -470,7 +417,9 @@ Foam::wallBoundedParticle::wallBoundedParticle const label diagEdge ) : - particle(mesh, position, celli, tetFacei, tetPti), +// particle(mesh, barycentric(1, 0, 0, 0), celli, tetFacei, tetPti), + particle(mesh, position, celli, tetFacei, tetPti, false), + localPosition_(position), meshEdgeStart_(meshEdgeStart), diagEdge_(diagEdge) {} @@ -490,11 +439,11 @@ Foam::wallBoundedParticle::wallBoundedParticle { if (is.format() == IOstream::ASCII) { - is >> meshEdgeStart_ >> diagEdge_; + is >> localPosition_ >> meshEdgeStart_ >> diagEdge_; } else { - is.read(reinterpret_cast(&meshEdgeStart_), sizeofFields_); + is.read(reinterpret_cast(&localPosition_), sizeofFields_); } } @@ -508,6 +457,7 @@ Foam::wallBoundedParticle::wallBoundedParticle ) : particle(p), + localPosition_(p.localPosition_), meshEdgeStart_(p.meshEdgeStart_), diagEdge_(p.diagEdge_) {} @@ -524,15 +474,16 @@ Foam::Ostream& Foam::operator<< if (os.format() == IOstream::ASCII) { os << static_cast(p) - << token::SPACE << p.meshEdgeStart_ - << token::SPACE << p.diagEdge_; + << token::SPACE << p.localPosition_ + << token::SPACE << p.meshEdgeStart_ + << token::SPACE << p.diagEdge_; } else { os << static_cast(p); os.write ( - reinterpret_cast(&p.meshEdgeStart_), + reinterpret_cast(&p.localPosition_), wallBoundedParticle::sizeofFields_ ); } @@ -568,7 +519,7 @@ Foam::Ostream& Foam::operator<< << " l 2 4" << nl << " l 3 4" << nl; os << " "; - meshTools::writeOBJ(os, p.position()); + meshTools::writeOBJ(os, p.localPosition_); return os; } diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H index c648fc5c3c..52f006ca1c 100644 --- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H +++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H @@ -71,10 +71,9 @@ class wallBoundedParticle public: //- Class used to pass tracking data to the trackToFace function - template - class TrackingData + class trackingData : - public particle::TrackingData + public particle::trackingData { public: @@ -83,16 +82,14 @@ public: // Constructors - inline TrackingData + template + inline trackingData ( - CloudType& cloud, + const TrackCloudType& cloud, const PackedBoolList& isWallPatch ) : - particle::TrackingData - ( - cloud - ), + particle::trackingData(cloud), isWallPatch_(isWallPatch) {} }; @@ -102,6 +99,10 @@ protected: // Protected data + //- Particle position is updated locally as opposed to via track + // functions of the base Foam::particle class + point localPosition_; + //- Particle is on mesh edge: // const face& f = mesh.faces()[tetFace()] // const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_)); @@ -123,10 +124,6 @@ protected: //- Construct current edge edge currentEdge() const; - //- Replacement for particle::currentTetIndices() that avoids bombing - // out on invalid tet decomposition (tetBasePtIs = -1) - tetIndices currentTetIndices() const; - //- Replacement for particle::crossEdgeConnectedFace that avoids bombing // out on invalid tet decomposition (tetBasePtIs = -1) void crossEdgeConnectedFace @@ -150,84 +147,6 @@ protected: bool isTriAlongTrack(const vector& n, const point& endPosition) const; - // Patch interactions - - //- Do all patch interaction - template - void patchInteraction(TrackData& td, const scalar trackFraction); - - //- Overridable function to handle the particle hitting a patch - // Executed before other patch-hitting functions - template - bool hitPatch - ( - const polyPatch&, - TrackData& td, - const label patchi, - const scalar trackFraction, - const tetIndices& tetIs - ); - - //- Overridable function to handle the particle hitting a wedge - template - void hitWedgePatch - ( - const wedgePolyPatch&, - TrackData& td - ); - - //- Overridable function to handle the particle hitting a - // symmetry plane - template - void hitSymmetryPlanePatch - ( - const symmetryPlanePolyPatch&, - TrackData& td - ); - - //- Overridable function to handle the particle hitting a - // symmetry patch - template - void hitSymmetryPatch - ( - const symmetryPolyPatch&, - TrackData& td - ); - - //- Overridable function to handle the particle hitting a cyclic - template - void hitCyclicPatch - ( - const cyclicPolyPatch&, - TrackData& td - ); - - //- Overridable function to handle the particle hitting a - //- processorPatch - template - void hitProcessorPatch - ( - const processorPolyPatch&, - TrackData& td - ); - - //- Overridable function to handle the particle hitting a wallPatch - template - void hitWallPatch - ( - const wallPolyPatch&, - TrackData& td, - const tetIndices& - ); - - //- Overridable function to handle the particle hitting a polyPatch - template - void hitPatch - ( - const polyPatch&, - TrackData& td - ); - public: // Constructors @@ -236,7 +155,7 @@ public: wallBoundedParticle ( const polyMesh& c, - const vector& position, + const point& position, const label celli, const label tetFacei, const label tetPti, @@ -308,14 +227,44 @@ public: // Track //- Equivalent of trackToFace - template + template scalar trackToEdge ( - TrackData& td, + TrackCloudType& cloud, + trackingData& td, const vector& endPosition ); + // Patch interactions + + //- Do all patch interaction + template + void patchInteraction + ( + TrackCloudType& cloud, + trackingData& td, + const scalar trackFraction + ); + + //- Overridable function to handle the particle hitting a + //- processorPatch + template + void hitProcessorPatch + ( + TrackCloudType& cloud, + trackingData& td + ); + + //- Overridable function to handle the particle hitting a wallPatch + template + void hitWallPatch + ( + TrackCloudType& cloud, + trackingData& td + ); + + // Info //- Return info proxy. diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C index b2febe5203..23f0a6044c 100644 --- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C +++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C @@ -27,97 +27,44 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -template +template void Foam::wallBoundedParticle::patchInteraction ( - TrackData& td, + TrackCloudType& cloud, + trackingData& td, const scalar trackFraction ) { - typedef typename TrackData::cloudType::particleType particleType; + typename TrackCloudType::particleType& p = + static_cast(*this); + typename TrackCloudType::particleType::trackingData& ttd = + static_cast(td); - particleType& p = static_cast(*this); - p.hitFace(td); - - if (!internalFace(facei_)) + if (!mesh().isInternalFace(face())) { - label origFacei = facei_; - label patchi = patch(facei_); + label origFacei = face(); + label patchi = patch(); - // 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 (face() != origFacei) { - // Did patch interaction model switch patches? - // Note: recalculate meshEdgeStart_, diagEdge_! - if (facei_ != origFacei) - { - patchi = patch(facei_); - } + patchi = patch(); + } - const polyPatch& patch = mesh_.boundaryMesh()[patchi]; + const polyPatch& patch = mesh().boundaryMesh()[patchi]; - if (isA(patch)) - { - p.hitWedgePatch - ( - static_cast(patch), td - ); - } - else if (isA(patch)) - { - p.hitSymmetryPlanePatch - ( - static_cast(patch), td - ); - } - else if (isA(patch)) - { - p.hitSymmetryPatch - ( - static_cast(patch), td - ); - } - else if (isA(patch)) - { - p.hitCyclicPatch - ( - static_cast(patch), td - ); - } - else if (isA(patch)) - { - p.hitProcessorPatch - ( - static_cast(patch), td - ); - } - else if (isA(patch)) - { - p.hitWallPatch - ( - static_cast(patch), td, faceHitTetIs - ); - } - else - { - p.hitPatch(patch, td); - } + if (isA(patch)) + { + p.hitProcessorPatch(cloud, ttd); + } + else if (isA(patch)) + { + p.hitWallPatch(cloud, ttd); + } + else + { + td.keepParticle = false; } } } @@ -125,10 +72,11 @@ void Foam::wallBoundedParticle::patchInteraction // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template +template Foam::scalar Foam::wallBoundedParticle::trackToEdge ( - TrackData& td, + TrackCloudType& cloud, + trackingData& td, const vector& endPosition ) { @@ -147,7 +95,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on) //checkInside(); - //checkOnTriangle(position()); + //checkOnTriangle(localPosition_); //if (meshEdgeStart_ != -1 || diagEdge_ != -1) //{ // checkOnEdge(); @@ -163,39 +111,36 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // 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())) + if (mesh().isInternalFace(tetFace())) { label nbrCelli = ( - celli_ == mesh_.faceOwner()[facei_] - ? mesh_.faceNeighbour()[facei_] - : mesh_.faceOwner()[facei_] + cell() == mesh().faceOwner()[face()] + ? mesh().faceNeighbour()[face()] + : mesh().faceOwner()[face()] ); // 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_); + tetIndices nbrTi(nbrCelli, tetFace(), tetPt()); - bool posVol = (nbrTi.tet(mesh_).mag() > 0); + const bool posVol = (nbrTi.tet(mesh()).mag() > 0); + const vector path(endPosition - localPosition_); - if - ( - posVol - == ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0) - ) + if (posVol == ((nbrTi.faceTri(mesh()).normal() & path) < 0)) { // Change into nbrCell. No need to change tetFace, tetPt. //Pout<< " crossed from cell:" << celli_ // << " into " << nbrCelli << endl; - celli_ = nbrCelli; - patchInteraction(td, trackFraction); + cell() = nbrCelli; + patchInteraction(cloud, td, trackFraction); } else { // Walk to other face on edge. Changes tetFace, tetPt but not // cell. crossEdgeConnectedFace(meshEdge); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } } else @@ -203,7 +148,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // Walk to other face on edge. This might give loop since // particle should have been removed? crossEdgeConnectedFace(meshEdge); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } } else @@ -211,41 +156,42 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // 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())) + if (mesh().isInternalFace(tetFace())) { FatalErrorInFunction << "Can only track on boundary faces." << " Face:" << tetFace() - << " at:" << mesh_.faceCentres()[tetFace()] + << " at:" << mesh().faceCentres()[tetFace()] << abort(FatalError); } - const triFace tri(currentTetIndices().faceTriIs(mesh_)); - vector n = tri.normal(mesh_.points()); + const triFace tri(currentTetIndices().faceTriIs(mesh(), false)); + vector n = tri.normal(mesh().points()); n /= mag(n); point projectedEndPosition = endPosition; - const bool posVol = (currentTetIndices().tet(mesh_).mag() > 0); + const bool posVol = (currentTetIndices().tet(mesh()).mag() > 0); if (!posVol) { // Negative tet volume. Track back by setting the end point - projectedEndPosition = position() - (endPosition-position()); + projectedEndPosition = + localPosition_ - (endPosition - localPosition_); // Make sure to use a large enough vector to cross the negative // face. Bit overkill. - const vector d(endPosition-position()); + const vector d(endPosition - localPosition_); const scalar magD(mag(d)); if (magD > ROOTVSMALL) { // Get overall mesh bounding box - treeBoundBox meshBb(mesh_.bounds()); + treeBoundBox meshBb(mesh().bounds()); // Extend to make 3D meshBb.inflate(ROOTSMALL); // Create vector guaranteed to cross mesh bounds - projectedEndPosition = position()-meshBb.mag()*d/magD; + projectedEndPosition = localPosition_ - meshBb.mag()*d/magD; // Clip to mesh bounds point intPt; @@ -253,9 +199,9 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge bool ok = meshBb.intersects ( projectedEndPosition, - position()-projectedEndPosition, + localPosition_ - projectedEndPosition, projectedEndPosition, - position(), + localPosition_, intPt, intPtBits ); @@ -269,8 +215,8 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // Remove normal component { - const point& basePt = mesh_.points()[tri[0]]; - projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n; + const point& basePt = mesh().points()[tri[0]]; + projectedEndPosition -= ((projectedEndPosition - basePt)&n)*n; } @@ -304,7 +250,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge } const tetIndices ti(currentTetIndices()); - + const triFace trif(ti.triIs(mesh(), false)); // Triangle (faceTriIs) gets constructed from // f[faceBasePtI_], // f[facePtAI_], @@ -315,21 +261,21 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // 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(); + const Foam::face& f = mesh().faces()[ti.face()]; + const label fp0 = trif[0]; if (triEdgei == 0) { - if (ti.facePtA() == f.fcIndex(fp0)) + if (trif[1] == f.fcIndex(fp0)) { //Pout<< "Real edge." << endl; diagEdge_ = -1; meshEdgeStart_ = fp0; //checkOnEdge(); crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } - else if (ti.facePtA() == f.rcIndex(fp0)) + else if (trif[1] == f.rcIndex(fp0)) { //Note: should not happen since boundary face so owner //Pout<< "Real edge." << endl; @@ -340,12 +286,12 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge meshEdgeStart_ = f.rcIndex(fp0); //checkOnEdge(); crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } else { // Get index of triangle on other side of edge. - diagEdge_ = ti.facePtA()-fp0; + diagEdge_ = trif[1] - fp0; if (diagEdge_ < 0) { diagEdge_ += f.size(); @@ -359,40 +305,39 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge { //Pout<< "Real edge." << endl; diagEdge_ = -1; - meshEdgeStart_ = ti.facePtA(); + meshEdgeStart_ = trif[1]; //checkOnEdge(); crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } else // if (triEdgei == 2) { - if (ti.facePtB() == f.rcIndex(fp0)) + if (trif[2] == f.rcIndex(fp0)) { //Pout<< "Real edge." << endl; diagEdge_ = -1; - meshEdgeStart_ = ti.facePtB(); + meshEdgeStart_ = trif[2]; //checkOnEdge(); crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } - else if (ti.facePtB() == f.fcIndex(fp0)) + else if (trif[2] == f.fcIndex(fp0)) { //Note: should not happen since boundary face so owner //Pout<< "Real edge." << endl; - FatalErrorInFunction - << abort(FatalError); + FatalErrorInFunction << abort(FatalError); diagEdge_ = -1; meshEdgeStart_ = fp0; //checkOnEdge(); crossEdgeConnectedFace(currentEdge()); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } else { //Pout<< "Triangle edge." << endl; // Get index of triangle on other side of edge. - diagEdge_ = ti.facePtB()-fp0; + diagEdge_ = trif[2] - fp0; if (diagEdge_ < 0) { diagEdge_ += f.size(); @@ -412,7 +357,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge // Particle is on mesh edge so change into other face on cell crossEdgeConnectedFace(currentEdge()); //checkOnEdge(); - patchInteraction(td, trackFraction); + patchInteraction(cloud, td, trackFraction); } else { @@ -430,74 +375,11 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge } -template -bool Foam::wallBoundedParticle::hitPatch -( - const polyPatch&, - TrackData& td, - const label patchi, - const scalar trackFraction, - const tetIndices& tetIs -) -{ - // Disable generic patch interaction - return false; -} - - -template -void Foam::wallBoundedParticle::hitWedgePatch -( - const wedgePolyPatch& pp, - TrackData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -template -void Foam::wallBoundedParticle::hitSymmetryPlanePatch -( - const symmetryPlanePolyPatch& pp, - TrackData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -template -void Foam::wallBoundedParticle::hitSymmetryPatch -( - const symmetryPolyPatch& pp, - TrackData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -template -void Foam::wallBoundedParticle::hitCyclicPatch -( - const cyclicPolyPatch& pp, - TrackData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -template +template void Foam::wallBoundedParticle::hitProcessorPatch ( - const processorPolyPatch& pp, - TrackData& td + TrackCloudType& cloud, + trackingData& td ) { // Switch particle @@ -508,44 +390,31 @@ void Foam::wallBoundedParticle::hitProcessorPatch // on the other side between 2 and 3 so edgeStart_ should be // f.size()-edgeStart_-1. - const Foam::face& f = mesh_.faces()[face()]; + const Foam::face& f = mesh().faces()[face()]; if (meshEdgeStart_ != -1) { - meshEdgeStart_ = f.size()-meshEdgeStart_-1; + meshEdgeStart_ = f.size() - meshEdgeStart_-1; } else { // diagEdge_ is relative to faceBasePt - diagEdge_ = f.size()-diagEdge_; + diagEdge_ = f.size() - diagEdge_; } } -template +template void Foam::wallBoundedParticle::hitWallPatch ( - const wallPolyPatch& wpp, - TrackData& td, - const tetIndices& + TrackCloudType& cloud, + trackingData& td ) {} -template -void Foam::wallBoundedParticle::hitPatch -( - const polyPatch& wpp, - TrackData& td -) -{ - // Remove particle - td.keepParticle = false; -} - - -template -void Foam::wallBoundedParticle::readFields(CloudType& c) +template +void Foam::wallBoundedParticle::readFields(TrackCloudType& c) { if (!c.size()) { @@ -554,34 +423,47 @@ void Foam::wallBoundedParticle::readFields(CloudType& c) particle::readFields(c); + IOField localPosition + ( + c.fieldIOobject("position", IOobject::MUST_READ) + ); + c.checkFieldIOobject(c, localPosition); + IOField