From e9fb8b8572750a358a405a72e99e4ad38863db43 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 1 Jun 2017 09:59:38 +0100 Subject: [PATCH] tetIndices: Removed duplicate logic The logic for generating tetrahedra from a face base point and an offset was duplicated in a few places. It is now confined to the tetIndices class. --- src/OpenFOAM/meshes/polyMesh/polyMesh.C | 11 +- .../polyMeshTetDecomposition.C | 129 +-------- .../polyMeshTetDecomposition.H | 9 - .../polyMeshTetDecomposition/tetIndices.C | 68 +---- .../polyMeshTetDecomposition/tetIndices.H | 119 +++------ .../polyMeshTetDecomposition/tetIndicesI.H | 244 +++++++----------- .../cfdTools/general/levelSet/levelSet.C | 159 ++++++++++++ .../general/levelSet/levelSetTemplates.C | 187 ++++++++++++++ .../cellPointWeight/cellPointWeight.C | 26 +- .../interpolationCellPointI.H | 15 +- src/lagrangian/basic/particle/particle.C | 91 ++----- src/lagrangian/basic/particle/particle.H | 22 -- src/lagrangian/basic/particle/particleI.H | 8 +- .../basic/particle/particleTemplates.C | 9 +- .../MPPICParcel/MPPICParcelTrackingDataI.H | 10 +- .../AveragingMethod/AveragingMethod.C | 13 +- .../MPPIC/AveragingMethods/Dual/Dual.C | 16 +- .../MPPIC/AveragingMethods/Moment/Moment.C | 11 +- .../DampingModels/Relaxation/Relaxation.C | 5 +- .../IsotropyModels/Stochastic/Stochastic.C | 10 +- .../MPPIC/PackingModels/Explicit/Explicit.C | 5 +- .../MPPIC/PackingModels/Implicit/Implicit.C | 4 +- 22 files changed, 560 insertions(+), 611 deletions(-) create mode 100644 src/finiteVolume/cfdTools/general/levelSet/levelSet.C create mode 100644 src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 52b8c4c54c..328fb417bd 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1389,16 +1389,7 @@ bool Foam::polyMesh::pointInCell for (label tetPti = 1; tetPti < f.size() - 1; tetPti++) { // Get tetIndices of face triangle - tetIndices faceTetIs - ( - polyMeshTetDecomposition::triangleTetIndices - ( - *this, - facei, - celli, - tetPti - ) - ); + tetIndices faceTetIs(celli, facei, tetPti); triPointRef faceTri = faceTetIs.faceTri(*this); diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C index 4f1edca365..7fa09a835d 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.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 @@ -530,11 +530,7 @@ Foam::List Foam::polyMeshTetDecomposition::faceTetIndices label cI ) { - static label nWarnings = 0; - static const label maxWarnings = 100; - const faceList& pFaces = mesh.faces(); - const labelList& pOwner = mesh.faceOwner(); const face& f = pFaces[fI]; @@ -542,123 +538,15 @@ Foam::List Foam::polyMeshTetDecomposition::faceTetIndices List faceTets(nTets); - bool own = (pOwner[fI] == cI); - - label tetBasePtI = mesh.tetBasePtIs()[fI]; - - if (tetBasePtI == -1) + for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++) { - if (nWarnings < maxWarnings) - { - WarningInFunction - << "No base point for face " << fI << ", " << f - << ", vertices " << UIndirectList(mesh.points(), f) - << ", produces a valid tet decomposition." - << endl; - nWarnings++; - } - if (nWarnings == maxWarnings) - { - Warning<< "Suppressing any further warnings." << endl; - nWarnings++; - } - - tetBasePtI = 0; - } - - for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++) - { - tetIndices& faceTetIs = faceTets[tetPtI - 1]; - - label facePtI = (tetPtI + tetBasePtI) % f.size(); - label otherFacePtI = f.fcIndex(facePtI); - - faceTetIs.cell() = cI; - - faceTetIs.face() = fI; - - faceTetIs.faceBasePt() = tetBasePtI; - - if (own) - { - faceTetIs.facePtA() = facePtI; - faceTetIs.facePtB() = otherFacePtI; - } - else - { - faceTetIs.facePtA() = otherFacePtI; - faceTetIs.facePtB() = facePtI; - } - - faceTetIs.tetPt() = tetPtI; + faceTets[tetPtI - 1] = tetIndices(cI, fI, tetPtI); } return faceTets; } -Foam::tetIndices Foam::polyMeshTetDecomposition::triangleTetIndices -( - const polyMesh& mesh, - const label fI, - const label cI, - const label tetPtI -) -{ - static label nWarnings = 0; - static const label maxWarnings = 100; - - const face& f = mesh.faces()[fI]; - bool own = (mesh.faceOwner()[fI] == cI); - label tetBasePtI = mesh.tetBasePtIs()[fI]; - if (tetBasePtI == -1) - { - if (nWarnings < maxWarnings) - { - WarningInFunction - << "No base point for face " << fI << ", " << f - << ", vertices " << UIndirectList(mesh.points(), f) - << ", produces a valid tet decomposition." - << endl; - nWarnings++; - } - if (nWarnings == maxWarnings) - { - Warning<< "Suppressing any further warnings." << endl; - nWarnings++; - } - - tetBasePtI = 0; - } - - tetIndices faceTetIs; - - label facePtI = (tetPtI + tetBasePtI) % f.size(); - label otherFacePtI = f.fcIndex(facePtI); - - faceTetIs.cell() = cI; - - faceTetIs.face() = fI; - - faceTetIs.faceBasePt() = tetBasePtI; - - if (own) - { - faceTetIs.facePtA() = facePtI; - faceTetIs.facePtB() = otherFacePtI; - } - else - { - faceTetIs.facePtA() = otherFacePtI; - faceTetIs.facePtB() = facePtI; - } - - faceTetIs.tetPt() = tetPtI; - - return faceTetIs; -} - - Foam::List Foam::polyMeshTetDecomposition::cellTetIndices ( const polyMesh& mesh, @@ -713,16 +601,7 @@ Foam::tetIndices Foam::polyMeshTetDecomposition::findTet for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++) { // Get tetIndices of face triangle - tetIndices faceTetIs - ( - triangleTetIndices - ( - mesh, - fI, - cI, - tetPtI - ) - ); + tetIndices faceTetIs(cI, fI, tetPtI); // Check if inside if (faceTetIs.tet(mesh).inside(pt)) diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H index 18de73777c..2a74b39da5 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H @@ -130,15 +130,6 @@ public: label cI ); - //- Return the tet decomposition of the given triangle of the given face - static tetIndices triangleTetIndices - ( - const polyMesh& mesh, - label fI, - label cI, - const label tetPtI // offset in face - ); - //- Return the tet decomposition of the given cell, see // findFacePt for the meaning of the indices static List cellTetIndices diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C index bfe4e2088e..f10c8ca695 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.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 @@ -25,15 +25,19 @@ License #include "tetIndices.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +const Foam::label Foam::tetIndices::maxNWarnings = 100; + +Foam::label Foam::tetIndices::nWarnings = 0; + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::tetIndices::tetIndices() : celli_(-1), facei_(-1), - faceBasePtI_(-1), - facePtAI_(-1), - facePtBI_(-1), tetPti_(-1) {} @@ -42,61 +46,15 @@ Foam::tetIndices::tetIndices ( label celli, label facei, - label faceBasePtI, - label facePtAI, - label facePtBI, label tetPtI ) : celli_(celli), facei_(facei), - faceBasePtI_(faceBasePtI), - facePtAI_(facePtAI), - facePtBI_(facePtBI), tetPti_(tetPtI) {} -Foam::tetIndices::tetIndices -( - label celli, - label facei, - label tetPtI, - const polyMesh& mesh -) -: - celli_(celli), - facei_(facei), - faceBasePtI_(-1), - facePtAI_(-1), - facePtBI_(-1), - tetPti_(tetPtI) -{ - const faceList& pFaces = mesh.faces(); - const labelList& pOwner = mesh.faceOwner(); - - const Foam::face& f = pFaces[facei_]; - - bool own = (pOwner[facei_] == celli_); - - faceBasePtI_ = mesh.tetBasePtIs()[facei_]; - - label facePtI = (tetPti_ + faceBasePtI_) % f.size(); - label otherFacePtI = f.fcIndex(facePtI); - - if (own) - { - facePtAI_ = facePtI; - facePtBI_ = otherFacePtI; - } - else - { - facePtAI_ = otherFacePtI; - facePtBI_ = facePtI; - } -} - - // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // Foam::tetIndices::~tetIndices() @@ -107,12 +65,7 @@ Foam::tetIndices::~tetIndices() Foam::Istream& Foam::operator>>(Istream& is, tetIndices& tI) { - is >> tI.cell() - >> tI.face() - >> tI.faceBasePt() - >> tI.facePtA() - >> tI.facePtB() - >> tI.tetPt(); + is >> tI.cell() >> tI.face() >> tI.tetPt(); is.check(FUNCTION_NAME); return is; @@ -123,9 +76,6 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const tetIndices& tI) { os << tI.cell() << token::SPACE << tI.face() << token::SPACE - << tI.faceBasePt() << token::SPACE - << tI.facePtA() << token::SPACE - << tI.facePtB() << token::SPACE << tI.tetPt() << token::SPACE << endl; diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H index 8ad320403a..8ab5d1e2e8 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H @@ -83,28 +83,25 @@ class tetIndices { // Private data - //- Cell that this is a decomposed tet of - label celli_; + //- Cell that this is a decomposed tet of + label celli_; - //- Face that holds this decomposed tet - label facei_; + //- Face that holds this decomposed tet + label facei_; - //- Base point on the face - label faceBasePtI_; + //- Point on the face, *relative to the base point*, which + // characterises this tet on the face. + label tetPti_; - //- Point on the face such that the right-hand circulation - // {faceBasePtI_, facePtAI_, facePtBI_} - // forms a triangle that points out of the tet - label facePtAI_; - //- Point on the face such that the right-hand circulation - // {faceBasePtI_, facePtAI_, facePtBI_} - // forms a triangle that points out of the tet - label facePtBI_; + // Private static data - //- Point on the face, *relative to the base point*, which - // characterises this tet on the face. - label tetPti_; + //- Maximum number of bad tet warnings + static const label maxNWarnings; + + //- Current number of bad tet warnings. Warnings stop when this reaches + // the maximum number. + static label nWarnings; public: @@ -115,24 +112,7 @@ public: tetIndices(); //- Construct from components - tetIndices - ( - label celli, - label facei, - label faceBasePtI, - label facePtAI, - label facePtBI, - label tetPtI - ); - - //- Construct from cell, face, pt description of tet - tetIndices - ( - label celli, - label facei, - label tetPtI, - const polyMesh& mesh - ); + tetIndices(label celli, label facei, label tetPtI); //- Destructor @@ -146,65 +126,36 @@ public: //- Return the cell inline label cell() const; - //- Return the face - inline label face() const; - - //- Return the face base point - inline label faceBasePt() const; - - //- Return face point A - inline label facePtA() const; - - //- Return face point B - inline label facePtB() const; - - //- Return the characterising tetPtI - inline label tetPt() const; - - //- Return the geometry corresponding to this tet from the - // supplied mesh - inline tetPointRef tet(const polyMesh& mesh) const; - - //- Return the geometry corresponding to this tet from the - // supplied mesh using the old positions - inline tetPointRef oldTet(const polyMesh& mesh) const; - - //- Return the geometry corresponding to the tri on the - // mesh face for this tet from the supplied mesh. Normal of - // the tri points out of the cell. - inline triPointRef faceTri(const polyMesh& mesh) const; - - //- Return the point indices corresponding to the tri on the mesh - // face for this tet from the supplied mesh. Normal of - // the tri points out of the cell. - inline triFace faceTriIs(const polyMesh& mesh) const; - - //- Return the geometry corresponding to the tri on the - // mesh face for this tet from the supplied mesh using - // the old position - inline triPointRef oldFaceTri(const polyMesh& mesh) const; - - - // Edit - //- Return non-const access to the cell inline label& cell(); + //- Return the face + inline label face() const; + //- Return non-const access to the face inline label& face(); - //- Return non-const access to the face base point - inline label& faceBasePt(); - - //- Return non-const access to face point A - inline label& facePtA(); - - //- Return non-const access to face point B - inline label& facePtB(); + //- Return the characterising tetPtI + inline label tetPt() const; //- Return non-const access to the characterising tetPtI 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; + + //- 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. + inline triPointRef faceTri(const polyMesh& mesh) const; + + //- Return the geometry corresponding to the tri on the face for + // this tet using the old positions. + inline triPointRef oldFaceTri(const polyMesh& mesh) const; + // Member Operators diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H index 1024595e4f..70edef09b1 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.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 @@ -25,181 +25,132 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::label Foam::tetIndices::cell() const +inline Foam::label Foam::tetIndices::cell() const { return celli_; } -Foam::label Foam::tetIndices::face() const -{ - return facei_; -} - - -Foam::label Foam::tetIndices::faceBasePt() const -{ - return faceBasePtI_; -} - - -Foam::label Foam::tetIndices::facePtA() const -{ - return facePtAI_; -} - - -Foam::label Foam::tetIndices::facePtB() const -{ - return facePtBI_; -} - - -Foam::label Foam::tetIndices::tetPt() const -{ - return tetPti_; -} - - -Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const -{ - const pointField& pPts = mesh.points(); - const faceList& pFaces = mesh.faces(); - const vectorField& pC = mesh.cellCentres(); - - const Foam::face& f = pFaces[facei_]; - - return tetPointRef - ( - pC[celli_], - pPts[f[faceBasePtI_]], - pPts[f[facePtAI_]], - pPts[f[facePtBI_]] - ); -} - - -Foam::tetPointRef Foam::tetIndices::oldTet(const polyMesh& mesh) const -{ - const pointField& oldPPts = mesh.oldPoints(); - const faceList& pFaces = mesh.faces(); - - // We need to reconstruct the old Cc from oldPoints (it isn't - // stored) - point oldC = mesh.cells()[celli_].centre - ( - oldPPts, - pFaces - ); - - const Foam::face& f = pFaces[facei_]; - - return tetPointRef - ( - oldC, - oldPPts[f[faceBasePtI_]], - oldPPts[f[facePtAI_]], - oldPPts[f[facePtBI_]] - ); -} - - -Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const -{ - const pointField& pPts = mesh.points(); - const faceList& pFaces = mesh.faces(); - - const Foam::face& f = pFaces[facei_]; - - return triPointRef - ( - pPts[f[faceBasePtI_]], - pPts[f[facePtAI_]], - pPts[f[facePtBI_]] - ); -} - - -Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const -{ - const faceList& pFaces = mesh.faces(); - - const Foam::face& f = pFaces[facei_]; - - return triFace - ( - f[faceBasePtI_], - f[facePtAI_], - f[facePtBI_] - ); -} - - -Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const -{ - const pointField& oldPPts = mesh.oldPoints(); - const faceList& pFaces = mesh.faces(); - - const Foam::face& f = pFaces[facei_]; - - return triPointRef - ( - oldPPts[f[faceBasePtI_]], - oldPPts[f[facePtAI_]], - oldPPts[f[facePtBI_]] - ); -} - - -Foam::label& Foam::tetIndices::cell() +inline Foam::label& Foam::tetIndices::cell() { return celli_; } -Foam::label& Foam::tetIndices::face() +inline Foam::label Foam::tetIndices::face() const { return facei_; } -Foam::label& Foam::tetIndices::faceBasePt() +inline Foam::label& Foam::tetIndices::face() { - return faceBasePtI_; + return facei_; } -Foam::label& Foam::tetIndices::facePtA() -{ - return facePtAI_; -} - - -Foam::label& Foam::tetIndices::facePtB() -{ - return facePtBI_; -} - - -Foam::label& Foam::tetIndices::tetPt() +inline Foam::label Foam::tetIndices::tetPt() const { return tetPti_; } +inline Foam::label& Foam::tetIndices::tetPt() +{ + return tetPti_; +} + + +inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const +{ + const Foam::face& f = mesh.faces()[face()]; + + label faceBasePtI = mesh.tetBasePtIs()[face()]; + + if (faceBasePtI < 0) + { + faceBasePtI = 0; + 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 faceOtherPtI = f.fcIndex(facePtI); + + if (mesh.faceOwner()[face()] != cell()) + { + Swap(facePtI, faceOtherPtI); + } + + return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]); +} + + +inline Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const +{ + const pointField& meshPoints = mesh.points(); + const triFace tri = faceTriIs(mesh); + + return tetPointRef + ( + mesh.cellCentres()[cell()], + meshPoints[tri[0]], + meshPoints[tri[1]], + meshPoints[tri[2]] + ); +} + + +inline Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const +{ + const pointField& meshPoints = mesh.points(); + const triFace tri = faceTriIs(mesh); + + return triPointRef + ( + meshPoints[tri[0]], + meshPoints[tri[1]], + meshPoints[tri[2]] + ); +} + + +inline Foam::triPointRef Foam::tetIndices::oldFaceTri +( + const polyMesh& mesh +) const +{ + const pointField& meshOldPoints = mesh.oldPoints(); + const triFace tri = faceTriIs(mesh); + + return triPointRef + ( + meshOldPoints[tri[0]], + meshOldPoints[tri[1]], + meshOldPoints[tri[2]] + ); +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // inline bool Foam::tetIndices::operator==(const Foam::tetIndices& rhs) const { return - ( cell() == rhs.cell() && face() == rhs.face() - && faceBasePt() == rhs.faceBasePt() - && facePtA() == rhs.facePtA() - && facePtB() == rhs.facePtB() - && tetPt() == rhs.tetPt() - ); + && tetPt() == rhs.tetPt(); } @@ -209,7 +160,4 @@ inline bool Foam::tetIndices::operator!=(const Foam::tetIndices& rhs) const } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - // ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C new file mode 100644 index 0000000000..e5a3fc9489 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C @@ -0,0 +1,159 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "levelSet.H" +#include "cut.H" +#include "polyMeshTetDecomposition.H" +#include "tetIndices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::tmp> +Foam::levelSetFraction +( + const fvMesh& mesh, + const scalarField& levelC, + const scalarField& levelP, + const bool above +) +{ + tmp> tResult + ( + new DimensionedField + ( + IOobject + ( + "levelSetFraction", + mesh.time().timeName(), + mesh + ), + mesh, + dimensionedScalar("0", dimless, 0) + ) + ); + DimensionedField& result = tResult.ref(); + + forAll(result, cI) + { + const List cellTetIs = + polyMeshTetDecomposition::cellTetIndices(mesh, cI); + + scalar v = 0, r = 0; + + forAll(cellTetIs, cellTetI) + { + const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh); + + const FixedList + tet = + { + mesh.cellCentres()[cI], + mesh.points()[triIs[0]], + mesh.points()[triIs[1]], + mesh.points()[triIs[2]] + }; + const FixedList + level = + { + levelC[cI], + levelP[triIs[0]], + levelP[triIs[1]], + levelP[triIs[2]] + }; + + v += cut::volumeOp()(tet); + + if (above) + { + r += tetCut(tet, level, cut::volumeOp(), cut::noOp()); + } + else + { + r += tetCut(tet, level, cut::noOp(), cut::volumeOp()); + } + } + + result[cI] = r/v; + } + + return tResult; +} + + +Foam::tmp Foam::levelSetFraction +( + const fvPatch& patch, + const scalarField& levelF, + const scalarField& levelP, + const bool above +) +{ + tmp tResult(new scalarField(patch.size(), 0)); + scalarField& result = tResult.ref(); + + forAll(result, fI) + { + const face& f = patch.patch().localFaces()[fI]; + + vector a = vector::zero, r = vector::zero; + + for(label eI = 0; eI < f.size(); ++ eI) + { + const edge e = f.faceEdge(eI); + + const FixedList + tri = + { + patch.patch().faceCentres()[fI], + patch.patch().localPoints()[e[0]], + patch.patch().localPoints()[e[1]] + }; + const FixedList + level = + { + levelF[fI], + levelP[e[0]], + levelP[e[1]] + }; + + a += cut::areaOp()(tri); + + if (above) + { + r += triCut(tri, level, cut::areaOp(), cut::noOp()); + } + else + { + r += triCut(tri, level, cut::noOp(), cut::areaOp()); + } + } + + result[fI] = a/magSqr(a) & r; + } + + return tResult; +} + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C new file mode 100644 index 0000000000..de89bbb76d --- /dev/null +++ b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C @@ -0,0 +1,187 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "levelSet.H" +#include "cut.H" +#include "polyMeshTetDecomposition.H" +#include "tetIndices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +Foam::tmp> Foam::levelSetAverage +( + const fvMesh& mesh, + const scalarField& levelC, + const scalarField& levelP, + const DimensionedField& positiveC, + const DimensionedField& positiveP, + const DimensionedField& negativeC, + const DimensionedField& negativeP +) +{ + tmp> tResult + ( + new DimensionedField + ( + IOobject + ( + positiveC.name() + ":levelSetAverage", + mesh.time().timeName(), + mesh + ), + mesh, + dimensioned("0", positiveC.dimensions(), Zero) + ) + ); + DimensionedField& result = tResult.ref(); + + forAll(result, cI) + { + const List cellTetIs = + polyMeshTetDecomposition::cellTetIndices(mesh, cI); + + scalar v = 0; + Type r = Zero; + + forAll(cellTetIs, cellTetI) + { + const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh); + + const FixedList + tet = + { + mesh.cellCentres()[cI], + mesh.points()[triIs[0]], + mesh.points()[triIs[1]], + mesh.points()[triIs[2]] + }; + const FixedList + level = + { + levelC[cI], + levelP[triIs[0]], + levelP[triIs[1]], + levelP[triIs[2]] + }; + const cut::volumeIntegrateOp + positive = FixedList + ({ + positiveC[cI], + positiveP[triIs[0]], + positiveP[triIs[1]], + positiveP[triIs[2]] + }); + const cut::volumeIntegrateOp + negative = FixedList + ({ + negativeC[cI], + negativeP[triIs[0]], + negativeP[triIs[1]], + negativeP[triIs[2]] + }); + + v += cut::volumeOp()(tet); + + r += tetCut(tet, level, positive, negative); + } + + result[cI] = r/v; + } + + return tResult; +} + + +template +Foam::tmp> Foam::levelSetAverage +( + const fvPatch& patch, + const scalarField& levelF, + const scalarField& levelP, + const Field& positiveF, + const Field& positiveP, + const Field& negativeF, + const Field& negativeP +) +{ + typedef typename outerProduct::type sumType; + + tmp> tResult(new Field(patch.size(), Zero)); + Field& result = tResult.ref(); + + forAll(result, fI) + { + const face& f = patch.patch().localFaces()[fI]; + + vector a = vector::zero; + sumType r = Zero; + + for(label eI = 0; eI < f.size(); ++ eI) + { + const edge e = f.faceEdge(eI); + + const FixedList + tri = + { + patch.patch().faceCentres()[fI], + patch.patch().localPoints()[e[0]], + patch.patch().localPoints()[e[1]] + }; + const FixedList + level = + { + levelF[fI], + levelP[e[0]], + levelP[e[1]] + }; + const cut::areaIntegrateOp + positive = FixedList + ({ + positiveF[fI], + positiveP[e[0]], + positiveP[e[1]] + }); + const cut::areaIntegrateOp + negative = FixedList + ({ + negativeF[fI], + negativeP[e[0]], + negativeP[e[1]] + }); + + a += cut::areaOp()(tri); + + r += triCut(tri, level, positive, negative); + } + + result[fI] = a/magSqr(a) & r; + } + + return tResult; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C index 218644c361..ab89551f84 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.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 @@ -55,15 +55,12 @@ void Foam::cellPointWeight::findTetrahedron celli ); - const faceList& pFaces = mesh.faces(); const scalar cellVolume = mesh.cellVolumes()[celli]; forAll(cellTets, tetI) { const tetIndices& tetIs = cellTets[tetI]; - const face& f = pFaces[tetIs.face()]; - // Barycentric coordinates of the position scalar det = tetIs.tet(mesh).barycentric(position, weights_); @@ -81,9 +78,8 @@ void Foam::cellPointWeight::findTetrahedron && (u + v + w < 1 + tol) ) { - faceVertices_[0] = f[tetIs.faceBasePt()]; - faceVertices_[1] = f[tetIs.facePtA()]; - faceVertices_[2] = f[tetIs.facePtB()]; + + faceVertices_ = tetIs.faceTriIs(mesh); return; } @@ -124,16 +120,12 @@ void Foam::cellPointWeight::findTetrahedron const tetIndices& tetIs = cellTets[nearestTetI]; - const face& f = pFaces[tetIs.face()]; - // Barycentric coordinates of the position, ignoring if the // determinant is suitable. If not, the return from barycentric // to weights_ is safe. tetIs.tet(mesh).barycentric(position, weights_); - faceVertices_[0] = f[tetIs.faceBasePt()]; - faceVertices_[1] = f[tetIs.facePtA()]; - faceVertices_[2] = f[tetIs.facePtB()]; + faceVertices_ = tetIs.faceTriIs(mesh); } @@ -160,8 +152,6 @@ void Foam::cellPointWeight::findTriangle const scalar faceAreaSqr = magSqr(mesh.faceAreas()[facei]); - const face& f = mesh.faces()[facei]; - forAll(faceTets, tetI) { const tetIndices& tetIs = faceTets[tetI]; @@ -189,9 +179,7 @@ void Foam::cellPointWeight::findTriangle weights_[2] = triWeights[1]; weights_[3] = triWeights[2]; - faceVertices_[0] = f[tetIs.faceBasePt()]; - faceVertices_[1] = f[tetIs.facePtA()]; - faceVertices_[2] = f[tetIs.facePtB()]; + faceVertices_ = tetIs.faceTriIs(mesh); return; } @@ -244,9 +232,7 @@ void Foam::cellPointWeight::findTriangle weights_[2] = triWeights[1]; weights_[3] = triWeights[2]; - faceVertices_[0] = f[tetIs.faceBasePt()]; - faceVertices_[1] = f[tetIs.facePtA()]; - faceVertices_[2] = f[tetIs.facePtB()]; + faceVertices_ = tetIs.faceTriIs(mesh); } diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H index 8cfa3947c8..e8e5bd3fa5 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.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 @@ -83,20 +83,15 @@ inline Type Foam::interpolationCellPoint::interpolate tetIs.tet(this->pMesh_).barycentric(position, weights); - const faceList& pFaces = this->pMesh_.faces(); - - const face& f = pFaces[tetIs.face()]; - // Order of weights is the same as that of the vertices of the tet, i.e. // cellCentre, faceBasePt, facePtA, facePtB. Type t = this->psi_[tetIs.cell()]*weights[0]; - t += psip_[f[tetIs.faceBasePt()]]*weights[1]; - - t += psip_[f[tetIs.facePtA()]]*weights[2]; - - t += psip_[f[tetIs.facePtB()]]*weights[3]; + const triFace triIs = tetIs.faceTriIs(this->pMesh_); + t += psip_[triIs[0]]*weights[1]; + t += psip_[triIs[1]]*weights[2]; + t += psip_[triIs[2]]*weights[3]; return t; } diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C index 6c62c85fbd..d206db410d 100644 --- a/src/lagrangian/basic/particle/particle.C +++ b/src/lagrangian/basic/particle/particle.C @@ -40,45 +40,6 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::particle::tetFaceIndices -( - label& baseI, - label& vertex1I, - label& vertex2I -) const -{ - const Foam::face& f = mesh_.faces()[tetFacei_]; - - baseI = max(0, mesh_.tetBasePtIs()[tetFacei_]); - - vertex1I = (baseI + tetPti_) % f.size(); - - vertex2I = f.fcIndex(vertex1I); - - if (mesh_.faceOwner()[tetFacei_] != celli_) - { - Swap(vertex1I, vertex2I); - } -} - - -void Foam::particle::tetMeshIndices -( - label& basei, - label& vertex1i, - label& vertex2i -) const -{ - const Foam::face& f = mesh_.faces()[tetFacei_]; - - tetFaceIndices(basei, vertex1i, vertex2i); - - basei = f[basei]; - vertex1i = f[vertex1i]; - vertex2i = f[vertex2i]; -} - - void Foam::particle::tetGeometry ( vector& centre, @@ -87,13 +48,14 @@ void Foam::particle::tetGeometry vector& vertex2 ) const { - label basei, vertex1i, vertex2i; - tetMeshIndices(basei, vertex1i, vertex2i); + const triFace triIs(currentTetIndices().faceTriIs(mesh_)); + const vectorField& ccs = mesh_.cellCentres(); + const pointField& pts = mesh_.points(); - centre = mesh_.cellCentres()[celli_]; - base = mesh_.points()[basei]; - vertex1 = mesh_.points()[vertex1i]; - vertex2 = mesh_.points()[vertex2i]; + centre = ccs[celli_]; + base = pts[triIs[0]]; + vertex1 = pts[triIs[1]]; + vertex2 = pts[triIs[2]]; } @@ -144,9 +106,7 @@ void Foam::particle::movingTetGeometry Pair& vertex2 ) const { - label basei, vertex1i, vertex2i; - tetMeshIndices(basei, vertex1i, vertex2i); - + const triFace triIs(currentTetIndices().faceTriIs(mesh_)); const pointField& ptsOld = mesh_.oldPoints(); const pointField& ptsNew = mesh_.points(); @@ -180,14 +140,14 @@ void Foam::particle::movingTetGeometry } centre[0] = ccOld + f0*(ccNew - ccOld); - base[0] = ptsOld[basei] + f0*(ptsNew[basei] - ptsOld[basei]); - vertex1[0] = ptsOld[vertex1i] + f0*(ptsNew[vertex1i] - ptsOld[vertex1i]); - vertex2[0] = ptsOld[vertex2i] + f0*(ptsNew[vertex2i] - ptsOld[vertex2i]); + 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[basei] - ptsOld[basei]); - vertex1[1] = f1*(ptsNew[vertex1i] - ptsOld[vertex1i]); - vertex2[1] = f1*(ptsNew[vertex2i] - ptsOld[vertex2i]); + 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]]); } @@ -364,23 +324,22 @@ void Foam::particle::changeTet(const label tetTriI) void Foam::particle::changeFace(const label tetTriI) { - // Get the tet topology - label basei, vertex1i, vertex2i; - tetMeshIndices(basei, vertex1i, vertex2i); + // Get the old topology + const triFace triOldIs(currentTetIndices().faceTriIs(mesh_)); // Get the shared edge and the pre-rotation edge sharedEdge; if (tetTriI == 1) { - sharedEdge = edge(vertex1i, vertex2i); + sharedEdge = edge(triOldIs[1], triOldIs[2]); } else if (tetTriI == 2) { - sharedEdge = edge(vertex2i, basei); + sharedEdge = edge(triOldIs[2], triOldIs[0]); } else if (tetTriI == 3) { - sharedEdge = edge(basei, vertex1i); + sharedEdge = edge(triOldIs[0], triOldIs[1]); } else { @@ -450,27 +409,27 @@ void Foam::particle::changeFace(const label tetTriI) } // Pre-rotation puts the shared edge opposite the base of the tetrahedron - if (sharedEdge.otherVertex(vertex1i) == -1) + if (sharedEdge.otherVertex(triOldIs[1]) == -1) { rotate(false); } - else if (sharedEdge.otherVertex(vertex2i) == -1) + else if (sharedEdge.otherVertex(triOldIs[2]) == -1) { rotate(true); } - // Update the new tet topology - tetMeshIndices(basei, vertex1i, vertex2i); + // Get the new topology + const triFace triNewIs = currentTetIndices().faceTriIs(mesh_); // Reflect to account for the change of triangle orientation on the new face reflect(); // Post rotation puts the shared edge back in the correct location - if (sharedEdge.otherVertex(vertex1i) == -1) + if (sharedEdge.otherVertex(triNewIs[1]) == -1) { rotate(true); } - else if (sharedEdge.otherVertex(vertex2i) == -1) + else if (sharedEdge.otherVertex(triNewIs[2]) == -1) { rotate(false); } diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 24355e3a59..341cf3de7f 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -174,24 +174,6 @@ private: // Tetrahedra functions - //- Get indices into the current face for the face-bound vertices of - // the current tet. - void tetFaceIndices - ( - label& baseI, - label& vertex1I, - label& vertex2I - ) const; - - //- Get indices into the mesh points for the face-bound vertices of - // the current tet. - void tetMeshIndices - ( - label& baseI, - label& vertex1I, - label& vertex2I - ) const; - //- Get the vertices of the current tet void tetGeometry ( @@ -482,10 +464,6 @@ public: // particle occupies. inline tetIndices currentTetIndices() const; - //- Return the geometry of the current tet that the - // particle occupies. - inline tetPointRef currentTet() const; - //- Return the normal of the tri on tetFacei_ for the // current tet. inline vector normal() const; diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H index 074bd8a5b0..fce5ded593 100644 --- a/src/lagrangian/basic/particle/particleI.H +++ b/src/lagrangian/basic/particle/particleI.H @@ -74,13 +74,7 @@ inline Foam::label Foam::particle::tetPt() const inline Foam::tetIndices Foam::particle::currentTetIndices() const { - return tetIndices(celli_, tetFacei_, tetPti_, mesh_); -} - - -inline Foam::tetPointRef Foam::particle::currentTet() const -{ - return currentTetIndices().tet(mesh_); + return tetIndices(celli_, tetFacei_, tetPti_); } diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index 748f093b3e..27291a03e8 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -218,14 +218,7 @@ void Foam::particle::trackToFace // No action is taken for tetPti_ for tetFacei_ here. These are handled // by the patch interaction call or later during processor transfer. - const tetIndices faceHitTetIs = - polyMeshTetDecomposition::triangleTetIndices - ( - mesh_, - tetFacei_, - celli_, - tetPti_ - ); + const tetIndices faceHitTetIs(celli_, tetFacei_, tetPti_); if ( diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H index 2be585ba41..1cebc01a23 100644 --- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H +++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -176,7 +176,7 @@ Foam::MPPICParcel::TrackingData::updateAverages forAllConstIter(typename CloudType, cloud, iter) { const typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh()); + const tetIndices tetIs = p.currentTetIndices(); const scalar m = p.nParticle()*p.mass(); @@ -194,7 +194,7 @@ Foam::MPPICParcel::TrackingData::updateAverages forAllConstIter(typename CloudType, cloud, iter) { const typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh()); + const tetIndices tetIs = p.currentTetIndices(); const vector u = uAverage_->interpolate(p.position(), tetIs); @@ -213,7 +213,7 @@ Foam::MPPICParcel::TrackingData::updateAverages forAllConstIter(typename CloudType, cloud, iter) { const typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh()); + const tetIndices tetIs = p.currentTetIndices(); weightAverage.add ( @@ -230,7 +230,7 @@ Foam::MPPICParcel::TrackingData::updateAverages forAllConstIter(typename CloudType, cloud, iter) { const typename CloudType::parcelType& p = iter(); - tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh()); + const tetIndices tetIs = p.currentTetIndices(); const scalar a = volumeAverage_->interpolate(p.position(), tetIs); const scalar r = radiusAverage_->interpolate(p.position(), tetIs); diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C index b26b326d6e..0a4820a814 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -203,20 +203,15 @@ bool Foam::AveragingMethod::write() const forAll(cellTets, tetI) { const tetIndices& tetIs = cellTets[tetI]; + const triFace triIs = tetIs.faceTriIs(mesh_); const scalar v = tetIs.tet(mesh_).mag(); cellValue[celli] += v*interpolate(mesh_.C()[celli], tetIs); cellGrad[celli] += v*interpolateGrad(mesh_.C()[celli], tetIs); - const face& f = mesh_.faces()[tetIs.face()]; - labelList vertices(3); - vertices[0] = f[tetIs.faceBasePt()]; - vertices[1] = f[tetIs.facePtA()]; - vertices[2] = f[tetIs.facePtB()]; - - forAll(vertices, vertexI) + forAll(triIs, vertexI) { - const label pointi = vertices[vertexI]; + const label pointi = triIs[vertexI]; pointVolume[pointi] += v; pointValue[pointi] += diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C index c0974c6a89..17fc6ba600 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -64,11 +64,11 @@ Foam::AveragingMethods::Dual::Dual forAll(cellTets, tetI) { const tetIndices& tetIs = cellTets[tetI]; - const face& f = this->mesh_.faces()[tetIs.face()]; + const triFace triIs = tetIs.faceTriIs(this->mesh_); const scalar v = tetIs.tet(this->mesh_).mag(); - volumeDual_[f[tetIs.faceBasePt()]] += v; - volumeDual_[f[tetIs.facePtA()]] += v; - volumeDual_[f[tetIs.facePtB()]] += v; + volumeDual_[triIs[0]] += v; + volumeDual_[triIs[1]] += v; + volumeDual_[triIs[2]] += v; } } @@ -113,11 +113,7 @@ void Foam::AveragingMethods::Dual::tetGeometry const tetIndices& tetIs ) const { - const face& f = this->mesh_.faces()[tetIs.face()]; - - tetVertices_[0] = f[tetIs.faceBasePt()]; - tetVertices_[1] = f[tetIs.facePtA()]; - tetVertices_[2] = f[tetIs.facePtB()]; + tetVertices_ = tetIs.faceTriIs(this->mesh_); tetIs.tet(this->mesh_).barycentric(position, tetCoordinates_); diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C index ba27159e9e..40f048d316 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,16 +70,15 @@ Foam::AveragingMethods::Moment::Moment forAll(cellTets, tetI) { const tetIndices& tetIs = cellTets[tetI]; - const label facei = tetIs.face(); - const face& f = mesh.faces()[facei]; + const triFace triIs = tetIs.faceTriIs(mesh); const tensor T ( tensor ( - mesh.points()[f[tetIs.faceBasePt()]] - mesh.C()[celli], - mesh.points()[f[tetIs.facePtA()]] - mesh.C()[celli], - mesh.points()[f[tetIs.facePtB()]] - mesh.C()[celli] + mesh.points()[triIs[0]] - mesh.C()[celli], + mesh.points()[triIs[1]] - mesh.C()[celli], + mesh.points()[triIs[2]] - mesh.C()[celli] ).T() ); diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C index d9f58b3b23..811c918edb 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -139,8 +139,7 @@ Foam::vector Foam::DampingModels::Relaxation::velocityCorrection const scalar deltaT ) const { - const tetIndices - tetIs(p.cell(), p.tetFace(), p.tetPt(), this->owner().mesh()); + const tetIndices tetIs(p.currentTetIndices()); const scalar x = deltaT*oneByTimeScaleAverage_->interpolate(p.position(), tetIs); diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C index 068568a1b5..cac1e4f9c3 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -166,7 +166,7 @@ void Foam::IsotropyModels::Stochastic::calculate() forAllIter(typename CloudType, this->owner(), iter) { typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); const scalar x = exponentAverage.interpolate(p.position(), tetIs); @@ -201,7 +201,7 @@ void Foam::IsotropyModels::Stochastic::calculate() forAllIter(typename CloudType, this->owner(), iter) { typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); uTildeAverage.add(p.position(), tetIs, p.nParticle()*p.mass()*p.U()); } uTildeAverage.average(massAverage); @@ -224,7 +224,7 @@ void Foam::IsotropyModels::Stochastic::calculate() forAllIter(typename CloudType, this->owner(), iter) { typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs); uTildeSqrAverage.add ( @@ -239,7 +239,7 @@ void Foam::IsotropyModels::Stochastic::calculate() forAllIter(typename CloudType, this->owner(), iter) { typename CloudType::parcelType& p = iter(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); const vector u = uAverage.interpolate(p.position(), tetIs); const scalar uRms = diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C index d11e55e35d..fda103663b 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -143,8 +143,7 @@ Foam::vector Foam::PackingModels::Explicit::velocityCorrection const scalar deltaT ) const { - const fvMesh& mesh = this->owner().mesh(); - const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); // interpolated quantities const scalar alpha = diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C index 33937408cd..9f82d53182 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -339,7 +339,7 @@ Foam::vector Foam::PackingModels::Implicit::velocityCorrection // containing tetrahedron and parcel coordinates within const label celli = p.cell(); const label facei = p.tetFace(); - const tetIndices tetIs(celli, facei, p.tetPt(), mesh); + const tetIndices tetIs(p.currentTetIndices()); FixedList tetCoordinates; tetIs.tet(mesh).barycentric(p.position(), tetCoordinates);