zhyang-dev/OpenFOAM-5.x
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

tetIndices: Removed duplicate logic

Browse Source 
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.
This commit is contained in:
Will Bainbridge
2017-06-01 09:59:38 +01:00
parent 1b377dd439
commit 3a85d0c91a
22 changed files with 234 additions and 639 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -1357,16 +1357,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);

127
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
View File

@ -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::tetIndices> 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,121 +538,15 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
List<tetIndices> 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
<< ", 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
<< ", 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::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
(
const polyMesh& mesh,
@ -711,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))

9
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
View File

@ -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<tetIndices> cellTetIndices

68
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C
View File

@ -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();
// Check state of Istream
is.check
@ -128,9 +81,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;

119
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View File

@ -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

244
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
View File

@ -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
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

19
src/finiteVolume/cfdTools/general/levelSet/levelSet.C
View File

@ -64,28 +64,23 @@ Foam::levelSetFraction
forAll(cellTetIs, cellTetI)
{
const tetIndices& tetIs = cellTetIs[cellTetI];
const face& f = mesh.faces()[tetIs.face()];
const label pI0 = f[tetIs.faceBasePt()];
const label pIA = f[tetIs.facePtA()];
const label pIB = f[tetIs.facePtB()];
const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
const FixedList<point, 4>
tet =
{
mesh.cellCentres()[cI],
mesh.points()[pI0],
mesh.points()[pIA],
mesh.points()[pIB]
mesh.points()[triIs[0]],
mesh.points()[triIs[1]],
mesh.points()[triIs[2]]
};
const FixedList<scalar, 4>
level =
{
levelC[cI],
levelP[pI0],
levelP[pIA],
levelP[pIB]
levelP[triIs[0]],
levelP[triIs[1]],
levelP[triIs[2]]
};
v += cut::volumeOp()(tet);

31
src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
View File

@ -68,44 +68,39 @@ Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>> Foam::levelSetAverage
forAll(cellTetIs, cellTetI)
{
const tetIndices& tetIs = cellTetIs[cellTetI];
const face& f = mesh.faces()[tetIs.face()];
const label pI0 = f[tetIs.faceBasePt()];
const label pIA = f[tetIs.facePtA()];
const label pIB = f[tetIs.facePtB()];
const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
const FixedList<point, 4>
tet =
{
mesh.cellCentres()[cI],
mesh.points()[pI0],
mesh.points()[pIA],
mesh.points()[pIB]
mesh.points()[triIs[0]],
mesh.points()[triIs[1]],
mesh.points()[triIs[2]]
};
const FixedList<scalar, 4>
level =
{
levelC[cI],
levelP[pI0],
levelP[pIA],
levelP[pIB]
levelP[triIs[0]],
levelP[triIs[1]],
levelP[triIs[2]]
};
const cut::volumeIntegrateOp<Type>
positive = FixedList<Type, 4>
({
positiveC[cI],
positiveP[pI0],
positiveP[pIA],
positiveP[pIB]
positiveP[triIs[0]],
positiveP[triIs[1]],
positiveP[triIs[2]]
});
const cut::volumeIntegrateOp<Type>
negative = FixedList<Type, 4>
({
negativeC[cI],
negativeP[pI0],
negativeP[pIA],
negativeP[pIB]
negativeP[triIs[0]],
negativeP[triIs[1]],
negativeP[triIs[2]]
});
v += cut::volumeOp()(tet);

26
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
View File

@ -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);
}

15
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H
View File

@ -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<Type>::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;
}

91
src/lagrangian/basic/particle/particle.C
View File

@ -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<vector>& 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);
}

22
src/lagrangian/basic/particle/particle.H
View File

@ -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
(
@ -476,10 +458,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;

8
src/lagrangian/basic/particle/particleI.H
View File

@ -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_);
}

9
src/lagrangian/basic/particle/particleTemplates.C
View File

@ -194,14 +194,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
(

10
src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
View File

@ -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<ParcelType>::TrackingData<CloudType>::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<ParcelType>::TrackingData<CloudType>::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<ParcelType>::TrackingData<CloudType>::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<ParcelType>::TrackingData<CloudType>::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);

13
src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
View File

@ -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
@ -204,20 +204,15 @@ bool Foam::AveragingMethod<Type>::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] +=

16
src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
View File

@ -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<Type>::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<Type>::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_);

11
src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
View File

@ -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<Type>::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()
);

5
src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
View File

@ -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<CloudType>::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);

10
src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
View File

@ -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<CloudType>::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<CloudType>::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<CloudType>::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<CloudType>::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 =

5
src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
View File

@ -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<CloudType>::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 =

4
src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
View File

@ -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<CloudType>::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<scalar, 4> tetCoordinates;
tetIs.tet(mesh).barycentric(p.position(), tetCoordinates);