From 81947c80859e6051aca73de4b766309fc7f0b60d Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 13 Mar 2018 08:22:33 +0000 Subject: [PATCH] particle: Optimisation by inlining key methods Resolves bug report https://bugs.openfoam.org/view.php?id=2871 --- src/lagrangian/basic/particle/particle.C | 124 ----------------- src/lagrangian/basic/particle/particle.H | 10 +- src/lagrangian/basic/particle/particleI.H | 128 +++++++++++++++++- .../basic/particle/particleTemplates.C | 9 +- 4 files changed, 135 insertions(+), 136 deletions(-) diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C index 508867f0bc..bc9dd3a840 100644 --- a/src/lagrangian/basic/particle/particle.C +++ b/src/lagrangian/basic/particle/particle.C @@ -42,34 +42,6 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::particle::stationaryTetGeometry -( - vector& centre, - vector& base, - vector& vertex1, - vector& vertex2 -) const -{ - const triFace triIs(currentTetIndices().faceTriIs(mesh_)); - const vectorField& ccs = mesh_.cellCentres(); - const pointField& pts = mesh_.points(); - - centre = ccs[celli_]; - base = pts[triIs[0]]; - vertex1 = pts[triIs[1]]; - vertex2 = pts[triIs[2]]; -} - - -Foam::barycentricTensor Foam::particle::stationaryTetTransform() const -{ - vector centre, base, vertex1, vertex2; - stationaryTetGeometry(centre, base, vertex1, vertex2); - - return barycentricTensor(centre, base, vertex1, vertex2); -} - - void Foam::particle::stationaryTetReverseTransform ( vector& centre, @@ -99,61 +71,6 @@ void Foam::particle::stationaryTetReverseTransform } -void Foam::particle::movingTetGeometry -( - const scalar fraction, - Pair& centre, - Pair& base, - Pair& vertex1, - Pair& vertex2 -) const -{ - const triFace triIs(currentTetIndices().faceTriIs(mesh_)); - const pointField& ptsOld = mesh_.oldPoints(); - const pointField& ptsNew = mesh_.points(); - - // !!! <-- We would be better off using mesh_.cellCentres() here. However, - // we need to put a mesh_.oldCellCentres() method in for this to work. The - // values obtained from the mesh and those obtained from the cell do not - // necessarily match. See mantis #1993. - const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces()); - const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces()); - - // 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. - 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]]); - vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]); - vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]); - - centre[1] = f1*(ccNew - ccOld); - base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]); - vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]); - vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]); -} - - -Foam::Pair Foam::particle::movingTetTransform -( - const scalar fraction -) const -{ - Pair centre, base, vertex1, vertex2; - movingTetGeometry(fraction, centre, base, vertex1, vertex2); - - return - Pair - ( - barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]), - barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1]) - ); -} - - void Foam::particle::movingTetReverseTransform ( const scalar fraction, @@ -961,47 +878,6 @@ Foam::vector Foam::particle::deviationFromMeshCentre() const } -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&) {} diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index c54f3c200c..ec89e957cf 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -168,7 +168,7 @@ private: // Tetrahedra functions //- Get the vertices of the current tet - void stationaryTetGeometry + inline void stationaryTetGeometry ( vector& centre, vector& base, @@ -181,7 +181,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 stationaryTetTransform() const; + inline 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 @@ -200,7 +200,7 @@ private: //- Get the vertices of the current moving tet. Two values are // returned for each vertex. The first is a constant, and the // second is a linear coefficient of the track fraction. - void movingTetGeometry + inline void movingTetGeometry ( const scalar endStepFraction, Pair& centre, @@ -213,7 +213,7 @@ private: // This is of the same form as for the static case. As with the // moving geometry, a linear function of the tracking fraction is // returned for each component. - Pair movingTetTransform + inline Pair movingTetTransform ( const scalar endStepFraction ) const; @@ -570,7 +570,7 @@ public: // Patch data //- Get the normal and velocity of the current patch location - void patchData(vector& n, vector& U) const; + inline 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 6d4c28a7cb..198b46c05a 100644 --- a/src/lagrangian/basic/particle/particleI.H +++ b/src/lagrangian/basic/particle/particleI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,6 +26,91 @@ License #include "polyMesh.H" #include "Time.H" +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::particle::stationaryTetGeometry +( + vector& centre, + vector& base, + vector& vertex1, + vector& vertex2 +) const +{ + const triFace triIs(currentTetIndices().faceTriIs(mesh_)); + const vectorField& ccs = mesh_.cellCentres(); + const pointField& pts = mesh_.points(); + + centre = ccs[celli_]; + base = pts[triIs[0]]; + vertex1 = pts[triIs[1]]; + vertex2 = pts[triIs[2]]; +} + + +inline Foam::barycentricTensor Foam::particle::stationaryTetTransform() const +{ + vector centre, base, vertex1, vertex2; + stationaryTetGeometry(centre, base, vertex1, vertex2); + + return barycentricTensor(centre, base, vertex1, vertex2); +} + + +inline void Foam::particle::movingTetGeometry +( + const scalar fraction, + Pair& centre, + Pair& base, + Pair& vertex1, + Pair& vertex2 +) const +{ + const triFace triIs(currentTetIndices().faceTriIs(mesh_)); + const pointField& ptsOld = mesh_.oldPoints(); + const pointField& ptsNew = mesh_.points(); + + // !!! <-- We would be better off using mesh_.cellCentres() here. However, + // we need to put a mesh_.oldCellCentres() method in for this to work. The + // values obtained from the mesh and those obtained from the cell do not + // necessarily match. See mantis #1993. + const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces()); + const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces()); + + // 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. + 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]]); + vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]); + vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]); + + centre[1] = f1*(ccNew - ccOld); + base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]); + vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]); + vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]); +} + + +inline Foam::Pair Foam::particle::movingTetTransform +( + const scalar fraction +) const +{ + Pair centre, base, vertex1, vertex2; + movingTetGeometry(fraction, centre, base, vertex1, vertex2); + + return + Pair + ( + barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]), + barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1]) + ); +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // inline Foam::label Foam::particle::getNewParticleID() const @@ -207,4 +292,45 @@ inline Foam::vector Foam::particle::position() const } +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; + } +} + + // ************************************************************************* // diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index 01f965c883..2b7ae6e012 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -122,7 +122,9 @@ void Foam::particle::hitFace } else if (onBoundaryFace()) { - if(!p.hitPatch(cloud, ttd)) + changeToMasterPatch(); + + if (!p.hitPatch(cloud, ttd)) { const polyPatch& patch = mesh_.boundaryMesh()[p.patch()]; @@ -178,11 +180,6 @@ void Foam::particle::trackToAndHitFace { trackToFace(direction, fraction); - if (onBoundaryFace()) - { - changeToMasterPatch(); - } - hitFace(direction, cloud, td); }