Lagrangian: Rewrite of the particle tracking algorithm to function in

terms of the local barycentric coordinates of the current tetrahedron,
rather than the global coordinate system.

Barycentric tracking works on any mesh, irrespective of mesh quality.
Particles do not get "lost", and tracking does not require ad-hoc
"corrections" or "rescues" to function robustly, because the calculation
of particle-face intersections is unambiguous and reproducible, even at
small angles of incidence.

Each particle position is defined by topology (i.e. the decomposed tet
cell it is in) and geometry (i.e. where it is in the cell). No search
operations are needed on restart or reconstruct, unlike when particle
positions are stored in the global coordinate system.

The particle positions file now contains particles' local coordinates
and topology, rather than the global coordinates and cell. This change
to the output format is not backwards compatible. Existing cases with
Lagrangian data will not restart, but they will still run from time
zero without any modification. This change was necessary in order to
guarantee that the loaded particle is valid, and therefore
fundamentally prevent "loss" and "search-failure" type bugs (e.g.,
2517, 2442, 2286, 1836, 1461, 1341, 1097).

The tracking functions have also been converted to function in terms
of displacement, rather than end position. This helps remove floating
point error issues, particularly towards the end of a tracking step.

Wall bounded streamlines have been removed. The implementation proved
incompatible with the new tracking algorithm. ParaView has a surface
LIC plugin which provides equivalent, or better, functionality.

Additionally, bug report <https://bugs.openfoam.org/view.php?id=2517>
is resolved by this change.
This commit is contained in:
Will Bainbridge
2017-04-28 08:03:44 +01:00
parent d2a62df7c4
commit 371762757d
99 changed files with 3081 additions and 5741 deletions

View File

@ -77,6 +77,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl; Info<< "Time = " << runTime.timeName() << nl << endl;
// Store the particle positions
kinematicCloud.storeGlobalPositions();
mesh.update(); mesh.update();
// Calculate absolute flux from the mapped surface velocity // Calculate absolute flux from the mapped surface velocity

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -65,6 +65,8 @@ int main(int argc, char *argv[])
{ {
Info<< "Time = " << runTime.timeName() << nl << endl; Info<< "Time = " << runTime.timeName() << nl << endl;
kinematicCloud.storeGlobalPositions();
mesh.update(); mesh.update();
U.correctBoundaryConditions(); U.correctBoundaryConditions();

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -88,6 +88,9 @@ int main(int argc, char *argv[])
// Store momentum to set rhoUf for introduced faces. // Store momentum to set rhoUf for introduced faces.
volVectorField rhoU("rhoU", rho*U); volVectorField rhoU("rhoU", rho*U);
// Store the particle positions
parcels.storeGlobalPositions();
// Do any mesh changes // Do any mesh changes
mesh.update(); mesh.update();

View File

@ -0,0 +1,3 @@
uncoupledKinematicParcelDyMFoam.C
EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelDyMFoam

View File

@ -0,0 +1,36 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-llagrangian \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lradiationModels \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-lregionModels \
-lsurfaceFilmModels \
-ldynamicMesh \
-ldynamicFvMesh \
-ltopoChangerFvMesh

View File

@ -0,0 +1,90 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Application
uncoupledKinematicParcelDyMFoam
Description
Transient solver for the passive transport of a particle cloud.
Uses a pre- calculated velocity field to evolve the cloud.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "psiThermo.H"
#include "turbulentFluidThermoModel.H"
#include "basicKinematicCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addOption
(
"cloudName",
"name",
"specify alternative cloud name. default is 'kinematicCloud'"
);
#define NO_CONTROL
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createFields.H"
#include "compressibleCourantNo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
kinematicCloud.storeGlobalPositions();
mesh.update();
U.correctBoundaryConditions();
Info<< "Evolving " << kinematicCloud.name() << endl;
kinematicCloud.evolve();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -47,13 +47,12 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
{ {
label pi = 0; label pi = 0;
// faceProcAddressing not required currently labelList decodedProcFaceAddressing(faceProcAddressing.size());
// labelList decodedProcFaceAddressing(faceProcAddressing.size());
// forAll(faceProcAddressing, i) forAll(faceProcAddressing, i)
// { {
// decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1; decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
// } }
forAll(cellProcAddressing, procCelli) forAll(cellProcAddressing, procCelli)
{ {
@ -68,27 +67,28 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
const indexedParticle& ppi = *iter(); const indexedParticle& ppi = *iter();
particleIndices_[pi++] = ppi.index(); particleIndices_[pi++] = ppi.index();
// label mappedTetFace = findIndex label mappedTetFace = findIndex
// ( (
// decodedProcFaceAddressing, decodedProcFaceAddressing,
// ppi.tetFace() ppi.tetFace()
// ); );
// if (mappedTetFace == -1) if (mappedTetFace == -1)
// { {
// FatalErrorInFunction FatalErrorInFunction
// << "Face lookup failure." << nl << "Face lookup failure." << nl
// << abort(FatalError); << abort(FatalError);
// } }
positions_.append positions_.append
( (
new passiveParticle new passiveParticle
( (
procMesh, procMesh,
ppi.position(), ppi.coordinates(),
procCelli, procCelli,
false mappedTetFace,
ppi.procTetPt(procMesh, procCelli, mappedTetFace)
) )
); );
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -98,7 +98,6 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
const fvMesh& meshSource = meshToMesh0Interp.fromMesh(); const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
const fvMesh& meshTarget = meshToMesh0Interp.toMesh(); const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs fileNameList cloudDirs
@ -182,15 +181,17 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
new passiveParticle new passiveParticle
( (
meshTarget, meshTarget,
targetCc[targetCells[i]], barycentric(1, 0, 0, 0),
targetCells[i] targetCells[i],
meshTarget.cells()[targetCells[i]][0],
1
) )
); );
passiveParticle& newP = newPtr(); passiveParticle& newP = newPtr();
label facei = newP.track(iter().position(), td); newP.track(iter().position() - newP.position(), 0);
if (facei < 0 && newP.cell() >= 0) if (!newP.onFace())
{ {
// Hit position. // Hit position.
foundCell = true; foundCell = true;
@ -233,11 +234,16 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
{ {
unmappedSource.erase(sourceParticleI); unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI); addParticles.append(sourceParticleI);
iter().cell() = targetCell;
targetParcels.addParticle targetParcels.addParticle
( (
sourceParcels.remove(&iter()) new passiveParticle
(
meshTarget,
iter().position(),
targetCell
)
); );
sourceParcels.remove(&iter());
} }
} }
sourceParticleI++; sourceParticleI++;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -89,8 +89,6 @@ void mapLagrangian(const meshToMesh& interp)
const polyMesh& meshTarget = interp.tgtRegion(); const polyMesh& meshTarget = interp.tgtRegion();
const labelListList& sourceToTarget = interp.srcToTgtCellAddr(); const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs fileNameList cloudDirs
( (
readDir readDir
@ -172,15 +170,17 @@ void mapLagrangian(const meshToMesh& interp)
new passiveParticle new passiveParticle
( (
meshTarget, meshTarget,
targetCc[targetCells[i]], barycentric(1, 0, 0, 0),
targetCells[i] targetCells[i],
meshTarget.cells()[targetCells[i]][0],
1
) )
); );
passiveParticle& newP = newPtr(); passiveParticle& newP = newPtr();
label facei = newP.track(iter().position(), td); newP.track(iter().position() - newP.position(), 0);
if (facei < 0 && newP.cell() >= 0) if (!newP.onFace())
{ {
// Hit position. // Hit position.
foundCell = true; foundCell = true;
@ -223,11 +223,16 @@ void mapLagrangian(const meshToMesh& interp)
{ {
unmappedSource.erase(sourceParticleI); unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI); addParticles.append(sourceParticleI);
iter().cell() = targetCell;
targetParcels.addParticle targetParcels.addParticle
( (
sourceParcels.remove(&iter()) new passiveParticle
(
meshTarget,
iter().position(),
targetCell
)
); );
sourceParcels.remove(&iter());
} }
} }
sourceParticleI++; sourceParticleI++;

View File

@ -33,8 +33,7 @@ minVol 1e-13;
//- Minimum quality of the tet formed by the face-centre //- Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and // and variable base point minimum decomposition triangles and
// the cell centre. This has to be a positive number for tracking // the cell centre. Set to very negative number (e.g. -1E30) to
// to work. Set to very negative number (e.g. -1E30) to
// disable. // disable.
// <0 = inside out tet, // <0 = inside out tet,
// 0 = flat tet // 0 = flat tet

View File

@ -128,6 +128,8 @@ $(spatialVectorAlgebra)/CompactSpatialTensor/compactSpatialTensor/compactSpatial
primitives/polynomialEqns/cubicEqn/cubicEqn.C primitives/polynomialEqns/cubicEqn/cubicEqn.C
primitives/polynomialEqns/quadraticEqn/quadraticEqn.C primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
primitives/Barycentric/barycentric/barycentric.C
containers/HashTables/HashTable/HashTableCore.C containers/HashTables/HashTable/HashTableCore.C
containers/HashTables/StaticHashTable/StaticHashTableCore.C containers/HashTables/StaticHashTable/StaticHashTableCore.C
containers/Lists/SortableList/ParSortableListName.C containers/Lists/SortableList/ParSortableListName.C

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -135,6 +135,29 @@ void Foam::polyMesh::calcDirections() const
} }
Foam::autoPtr<Foam::labelIOList> Foam::polyMesh::readTetBasePtIs() const
{
IOobject io
(
"tetBasePtIs",
instance(),
meshSubDir,
*this,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
);
if (io.headerOk())
{
return autoPtr<labelIOList>(new labelIOList(io));
}
else
{
return autoPtr<labelIOList>(nullptr);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyMesh::polyMesh(const IOobject& io) Foam::polyMesh::polyMesh(const IOobject& io)
@ -207,7 +230,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
comm_(UPstream::worldComm), comm_(UPstream::worldComm),
geometricD_(Zero), geometricD_(Zero),
solutionD_(Zero), solutionD_(Zero),
tetBasePtIsPtr_(nullptr), tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr), cellTreePtr_(nullptr),
pointZones_ pointZones_
( (
@ -401,7 +424,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm), comm_(UPstream::worldComm),
geometricD_(Zero), geometricD_(Zero),
solutionD_(Zero), solutionD_(Zero),
tetBasePtIsPtr_(nullptr), tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr), cellTreePtr_(nullptr),
pointZones_ pointZones_
( (
@ -552,7 +575,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm), comm_(UPstream::worldComm),
geometricD_(Zero), geometricD_(Zero),
solutionD_(Zero), solutionD_(Zero),
tetBasePtIsPtr_(nullptr), tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr), cellTreePtr_(nullptr),
pointZones_ pointZones_
( (
@ -815,7 +838,7 @@ Foam::label Foam::polyMesh::nSolutionD() const
} }
const Foam::labelList& Foam::polyMesh::tetBasePtIs() const const Foam::labelIOList& Foam::polyMesh::tetBasePtIs() const
{ {
if (tetBasePtIsPtr_.empty()) if (tetBasePtIsPtr_.empty())
{ {
@ -828,8 +851,17 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
tetBasePtIsPtr_.reset tetBasePtIsPtr_.reset
( (
new labelList new labelIOList
( (
IOobject
(
"tetBasePtIs",
instance(),
meshSubDir,
*this,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
polyMeshTetDecomposition::findFaceBasePts(*this) polyMeshTetDecomposition::findFaceBasePts(*this)
) )
); );
@ -1088,6 +1120,13 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
points_.instance() = time().timeName(); points_.instance() = time().timeName();
points_.eventNo() = getEvent(); points_.eventNo() = getEvent();
if (tetBasePtIsPtr_.valid())
{
tetBasePtIsPtr_().writeOpt() = IOobject::AUTO_WRITE;
tetBasePtIsPtr_().instance() = time().timeName();
tetBasePtIsPtr_().eventNo() = getEvent();
}
tmp<scalarField> sweptVols = primitiveMesh::movePoints tmp<scalarField> sweptVols = primitiveMesh::movePoints
( (
points_, points_,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -150,7 +150,7 @@ private:
mutable Vector<label> solutionD_; mutable Vector<label> solutionD_;
//- Base point for face decomposition into tets //- Base point for face decomposition into tets
mutable autoPtr<labelList> tetBasePtIsPtr_; mutable autoPtr<labelIOList> tetBasePtIsPtr_;
//- Search tree to allow spatial cell searching //- Search tree to allow spatial cell searching
mutable autoPtr<indexedOctree<treeDataCell>> cellTreePtr_; mutable autoPtr<indexedOctree<treeDataCell>> cellTreePtr_;
@ -208,6 +208,9 @@ private:
// polyhedral information // polyhedral information
void calcCellShapes() const; void calcCellShapes() const;
//- Read and return the tetBasePtIs
autoPtr<labelIOList> readTetBasePtIs() const;
// Helper functions for constructor from cell shapes // Helper functions for constructor from cell shapes
@ -449,7 +452,7 @@ public:
label nSolutionD() const; label nSolutionD() const;
//- Return the tetBasePtIs //- Return the tetBasePtIs
const labelList& tetBasePtIs() const; const labelIOList& tetBasePtIs() const;
//- Return the cell search tree //- Return the cell search tree
const indexedOctree<treeDataCell>& cellTree() const; const indexedOctree<treeDataCell>& cellTree() const;
@ -605,8 +608,8 @@ public:
//- Clear primitive data (points, faces and cells) //- Clear primitive data (points, faces and cells)
void clearPrimitives(); void clearPrimitives();
//- Clear geometry not used for CFD (cellTree, tetBasePtIs) //- Clear tet base points
void clearAdditionalGeom(); void clearTetBasePtIs();
//- Clear cell tree data //- Clear cell tree data
void clearCellTree(); void clearCellTree();

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -69,22 +69,6 @@ void Foam::polyMesh::clearGeom()
geometricD_ = Zero; geometricD_ = Zero;
solutionD_ = Zero; solutionD_ = Zero;
// Remove the stored tet base points
tetBasePtIsPtr_.clear();
// Remove the cell tree
cellTreePtr_.clear();
}
void Foam::polyMesh::clearAdditionalGeom()
{
if (debug)
{
InfoInFunction << "Clearing additional geometric data" << endl;
}
// Remove the stored tet base points
tetBasePtIsPtr_.clear();
// Remove the cell tree // Remove the cell tree
cellTreePtr_.clear(); cellTreePtr_.clear();
} }
@ -170,6 +154,17 @@ void Foam::polyMesh::clearOut()
} }
void Foam::polyMesh::clearTetBasePtIs()
{
if (debug)
{
InfoInFunction << "Clearing tet base points" << endl;
}
tetBasePtIsPtr_.clear();
}
void Foam::polyMesh::clearCellTree() void Foam::polyMesh::clearCellTree()
{ {
if (debug) if (debug)

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -59,6 +59,12 @@ void Foam::polyMesh::setInstance(const fileName& inst)
cellZones_.writeOpt() = IOobject::AUTO_WRITE; cellZones_.writeOpt() = IOobject::AUTO_WRITE;
cellZones_.instance() = inst; cellZones_.instance() = inst;
if (tetBasePtIsPtr_.valid())
{
tetBasePtIsPtr_->writeOpt() = IOobject::AUTO_WRITE;
tetBasePtIsPtr_->instance() = inst;
}
} }
@ -386,6 +392,9 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
cellZones_.set(czI, newCellZones[czI].clone(cellZones_)); cellZones_.set(czI, newCellZones[czI].clone(cellZones_));
} }
// Re-read tet base points
tetBasePtIsPtr_ = readTetBasePtIs();
if (boundaryChanged) if (boundaryChanged)
{ {
@ -438,6 +447,13 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
points_.transfer(newPoints); points_.transfer(newPoints);
points_.instance() = pointsInst; points_.instance() = pointsInst;
// Re-read tet base points
autoPtr<labelIOList> newTetBasePtIsPtr = readTetBasePtIs();
if (newTetBasePtIsPtr.valid())
{
tetBasePtIsPtr_ = newTetBasePtIsPtr;
}
// Calculate the geometry for the patches (transformation tensors etc.) // Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry(); boundary_.calcGeometry();

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -22,21 +22,21 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class Class
Foam::wallBoundedStreamLineParticleCloud Foam::Barycentric
Description Description
A Cloud of streamLine particles Templated 3D Barycentric derived from VectorSpace. Has 4 components, one of
which is redundant.
SourceFiles SourceFiles
streamLineCloud.C BarycentricI.H
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef streamLineCloud_H #ifndef Barycentric_H
#define streamLineCloud_H #define Barycentric_H
#include "Cloud.H" #include "VectorSpace.H"
#include "wallBoundedStreamLineParticle.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,47 +44,61 @@ namespace Foam
{ {
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class streamLineCloud Declaration Class Barycentric Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
class wallBoundedStreamLineParticleCloud template<class Cmpt>
class Barycentric
: :
public Cloud<wallBoundedStreamLineParticle> public VectorSpace<Barycentric<Cmpt>, Cmpt, 4>
{ {
// Private Member Functions
//- Disallow default bitwise copy construct
wallBoundedStreamLineParticleCloud
(
const wallBoundedStreamLineParticleCloud&
);
//- Disallow default bitwise assignment
void operator=(const wallBoundedStreamLineParticleCloud&);
public: public:
//- Type of parcel the cloud was instantiated for //- Equivalent type of labels used for valid component indexing
typedef wallBoundedStreamLineParticle parcelType; typedef Barycentric<label> labelType;
// Member constants
//- Rank of Barycentric is 1
static const direction rank = 1;
//- Component labeling enumeration
enum components { A, B, C, D };
// Constructors // Constructors
//- Construct given mesh //- Construct null
wallBoundedStreamLineParticleCloud inline Barycentric();
//- Construct initialized to zero
inline Barycentric(const Foam::zero);
//- Construct given four components
inline Barycentric
( (
const polyMesh&, const Cmpt& va,
const word& cloudName = "defaultCloud", const Cmpt& vb,
bool readFields = true const Cmpt& vc,
const Cmpt& vd
); );
//- Construct from mesh, cloud name, and a list of particles
wallBoundedStreamLineParticleCloud // Member Functions
(
const polyMesh& mesh, // Access
const word& cloudName,
const IDLList<wallBoundedStreamLineParticle>& particles inline const Cmpt& a() const;
); inline const Cmpt& b() const;
inline const Cmpt& c() const;
inline const Cmpt& d() const;
inline Cmpt& a();
inline Cmpt& b();
inline Cmpt& c();
inline Cmpt& d();
}; };
@ -94,6 +108,10 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "BarycentricI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif #endif
// ************************************************************************* // // ************************************************************************* //

View File

@ -0,0 +1,132 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-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
ANB 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::Barycentric<Cmpt>::Barycentric()
{}
template<class Cmpt>
inline Foam::Barycentric<Cmpt>::Barycentric(const Foam::zero)
:
Barycentric::vsType(Zero)
{}
template<class Cmpt>
inline Foam::Barycentric<Cmpt>::Barycentric
(
const Cmpt& va,
const Cmpt& vb,
const Cmpt& vc,
const Cmpt& vd
)
{
this->v_[A] = va;
this->v_[B] = vb;
this->v_[C] = vc;
this->v_[D] = vd;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline const Cmpt& Foam::Barycentric<Cmpt>::a() const
{
return this->v_[A];
}
template<class Cmpt>
inline const Cmpt& Foam::Barycentric<Cmpt>::b() const
{
return this->v_[B];
}
template<class Cmpt>
inline const Cmpt& Foam::Barycentric<Cmpt>::c() const
{
return this->v_[C];
}
template<class Cmpt>
inline const Cmpt& Foam::Barycentric<Cmpt>::d() const
{
return this->v_[D];
}
template<class Cmpt>
inline Cmpt& Foam::Barycentric<Cmpt>::a()
{
return this->v_[A];
}
template<class Cmpt>
inline Cmpt& Foam::Barycentric<Cmpt>::b()
{
return this->v_[B];
}
template<class Cmpt>
inline Cmpt& Foam::Barycentric<Cmpt>::c()
{
return this->v_[C];
}
template<class Cmpt>
inline Cmpt& Foam::Barycentric<Cmpt>::d()
{
return this->v_[D];
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class Cmpt>
inline Cmpt operator&(const Barycentric<Cmpt>& b1, const Barycentric<Cmpt>& b2)
{
return b1.a()*b2.a() + b1.b()*b2.b() + b1.c()*b2.c() + b1.d()*b2.d();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,139 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
Foam::BarycentricTensor
Description
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can
represent a barycentric transformation as a matrix-barycentric inner-
product. Can alternatively represent an inverse barycentric transformation
as a vector-matrix inner-product.
SourceFiles
BarycentricTensorI.H
\*---------------------------------------------------------------------------*/
#ifndef BarycentricTensor_H
#define BarycentricTensor_H
#include "Barycentric.H"
#include "Tensor.H"
#include "Vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class BarycentricTensor Declaration
\*---------------------------------------------------------------------------*/
template<class Cmpt>
class BarycentricTensor
:
public MatrixSpace<BarycentricTensor<Cmpt>, Cmpt, 4, 3>
{
public:
//- Equivalent type of labels used for valid component indexing
typedef Tensor<label> labelType;
// Member constants
//- Rank of BarycentricTensor is 2
static const direction rank = 2;
//- Component labeling enumeration
enum components { XA, XB, XC, XD, YA, YB, YC, YD, ZA, ZB, ZC, ZD };
// Constructors
//- Construct null
BarycentricTensor();
//- Construct initialised to zero
BarycentricTensor(const Foam::zero);
//- Construct given three barycentric components (rows)
BarycentricTensor
(
const Barycentric<Cmpt>& x,
const Barycentric<Cmpt>& y,
const Barycentric<Cmpt>& z
);
//- Construct given four vector components (columns)
BarycentricTensor
(
const Vector<Cmpt>& a,
const Vector<Cmpt>& b,
const Vector<Cmpt>& c,
const Vector<Cmpt>& d
);
// Member Functions
// Row-barycentric access
inline Barycentric<Cmpt> x() const;
inline Barycentric<Cmpt> y() const;
inline Barycentric<Cmpt> z() const;
// Column-vector access
inline Vector<Cmpt> a() const;
inline Vector<Cmpt> b() const;
inline Vector<Cmpt> c() const;
inline Vector<Cmpt> d() const;
};
template<class Cmpt>
class typeOfTranspose<Cmpt, BarycentricTensor<Cmpt>>
{
public:
typedef void type;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "BarycentricTensorI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,196 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor()
{}
template<class Cmpt>
inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor(const Foam::zero)
:
BarycentricTensor::msType(Zero)
{}
template<class Cmpt>
inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor
(
const Barycentric<Cmpt>& x,
const Barycentric<Cmpt>& y,
const Barycentric<Cmpt>& z
)
{
this->v_[XA] = x.a();
this->v_[XB] = x.b();
this->v_[XC] = x.c();
this->v_[XD] = x.d();
this->v_[YA] = y.a();
this->v_[YB] = y.b();
this->v_[YC] = y.c();
this->v_[YD] = y.d();
this->v_[ZA] = z.a();
this->v_[ZB] = z.b();
this->v_[ZC] = z.c();
this->v_[ZD] = z.d();
}
template<class Cmpt>
inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor
(
const Vector<Cmpt>& a,
const Vector<Cmpt>& b,
const Vector<Cmpt>& c,
const Vector<Cmpt>& d
)
{
this->v_[XA] = a.x();
this->v_[XB] = b.x();
this->v_[XC] = c.x();
this->v_[XD] = d.x();
this->v_[YA] = a.y();
this->v_[YB] = b.y();
this->v_[YC] = c.y();
this->v_[YD] = d.y();
this->v_[ZA] = a.z();
this->v_[ZB] = b.z();
this->v_[ZC] = c.z();
this->v_[ZD] = d.z();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::x() const
{
return
Barycentric<Cmpt>
(
this->v_[XA],
this->v_[XB],
this->v_[XC],
this->v_[XD]
);
}
template<class Cmpt>
inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::y() const
{
return
Barycentric<Cmpt>
(
this->v_[YA],
this->v_[YB],
this->v_[YC],
this->v_[YD]
);
}
template<class Cmpt>
inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::z() const
{
return
Barycentric<Cmpt>
(
this->v_[ZA],
this->v_[ZB],
this->v_[ZC],
this->v_[ZD]
);
}
template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::a() const
{
return Vector<Cmpt>(this->v_[XA], this->v_[YA], this->v_[ZA]);
}
template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::b() const
{
return Vector<Cmpt>(this->v_[XB], this->v_[YB], this->v_[ZB]);
}
template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::c() const
{
return Vector<Cmpt>(this->v_[XC], this->v_[YC], this->v_[ZC]);
}
template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::d() const
{
return Vector<Cmpt>(this->v_[XD], this->v_[YD], this->v_[ZD]);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class Cmpt>
inline Vector<Cmpt> operator&
(
const BarycentricTensor<Cmpt>& T,
const Barycentric<Cmpt>& b
)
{
return Vector<Cmpt>(T.x() & b, T.y() & b, T.z() & b);
}
template<class Cmpt>
inline Barycentric<Cmpt> operator&
(
const Vector<Cmpt>& v,
const BarycentricTensor<Cmpt>& T
)
{
return Barycentric<Cmpt>(v & T.a(), v & T.b(), v & T.c(), v & T.d());
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,95 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "barycentric.H"
#include "Random.H"
#include "cachedRandom.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::barycentric barycentric01
(
Foam::scalar s,
Foam::scalar t,
Foam::scalar u
)
{
// Transform the random point in the unit cube to a random point in the
// unit tet by means of a series of reflections. See
// <http://vcg.isti.cnr.it/jgt/tetra.htm> for details.
if (s + t > 1)
{
s = 1 - s;
t = 1 - t;
}
if (s + t + u > 1)
{
Foam::scalar temp = u;
if (t + u > 1)
{
u = 1 - s - t;
t = 1 - temp;
}
else
{
u = s + t + u - 1;
s = 1 - t - temp;
}
}
return Foam::barycentric(1 - s - t - u, s, t, u);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::barycentric Foam::barycentric01(Random& rndGen)
{
return
::barycentric01
(
rndGen.scalar01(),
rndGen.scalar01(),
rndGen.scalar01()
);
}
Foam::barycentric Foam::barycentric01(cachedRandom& rndGen)
{
return
::barycentric01
(
rndGen.sample01<scalar>(),
rndGen.sample01<scalar>(),
rndGen.sample01<scalar>()
);
}
// ************************************************************************* //

View File

@ -0,0 +1,74 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-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 <http://www.gnu.org/licenses/>.
Typedef
Foam::barycentric
Description
A scalar version of the templated Barycentric
\*---------------------------------------------------------------------------*/
#ifndef barycentric_H
#define barycentric_H
#include "scalar.H"
#include "Barycentric.H"
#include "contiguous.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Random;
class cachedRandom;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
typedef Barycentric<scalar> barycentric;
//- Generate a random barycentric coordinate within the unit tetrahedron
barycentric barycentric01(Random& rndGen);
barycentric barycentric01(cachedRandom& rndGen);
template<>
inline bool contiguous<barycentric>()
{
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -21,44 +21,44 @@ License
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::barycentricTensor
Description
A scalar version of the templated BarycentricTensor
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "wallBoundedStreamLineParticleCloud.H" #ifndef barycentricTensor_H
#define barycentricTensor_H
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // #include "scalar.H"
#include "BarycentricTensor.H"
#include "contiguous.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam namespace Foam
{ {
defineTemplateTypeNameAndDebug(Cloud<wallBoundedStreamLineParticle>, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud typedef BarycentricTensor<scalar> barycentricTensor;
(
const polyMesh& mesh,
const word& cloudName, template<>
bool readFields inline bool contiguous<barycentricTensor>()
)
:
Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, false)
{ {
if (readFields) return true;
{
wallBoundedStreamLineParticle::readFields(*this);
}
} }
Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
(
const polyMesh& mesh,
const word& cloudName,
const IDLList<wallBoundedStreamLineParticle>& particles
)
:
Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, particles)
{}
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* // // ************************************************************************* //

View File

@ -92,6 +92,9 @@ public:
//- Evaluate at x //- Evaluate at x
inline scalar value(const scalar x) const; inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x //- Estimate the error of evaluation at x
inline scalar error(const scalar x) const; inline scalar error(const scalar x) const;

View File

@ -106,6 +106,12 @@ inline Foam::scalar Foam::cubicEqn::value(const scalar x) const
} }
inline Foam::scalar Foam::cubicEqn::derivative(const scalar x) const
{
return x*(x*3*a() + 2*b()) + c();
}
inline Foam::scalar Foam::cubicEqn::error(const scalar x) const inline Foam::scalar Foam::cubicEqn::error(const scalar x) const
{ {
return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c())); return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c()));

View File

@ -81,6 +81,9 @@ public:
//- Evaluate at x //- Evaluate at x
inline scalar value(const scalar x) const; inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x //- Estimate the error of evaluation at x
inline scalar error(const scalar x) const; inline scalar error(const scalar x) const;

View File

@ -74,6 +74,12 @@ inline Foam::scalar Foam::linearEqn::value(const scalar x) const
} }
inline Foam::scalar Foam::linearEqn::derivative(const scalar x) const
{
return a();
}
inline Foam::scalar Foam::linearEqn::error(const scalar x) const inline Foam::scalar Foam::linearEqn::error(const scalar x) const
{ {
return mag(SMALL*x*a()); return mag(SMALL*x*a());

View File

@ -84,6 +84,9 @@ public:
//- Evaluate at x //- Evaluate at x
inline scalar value(const scalar x) const; inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x //- Estimate the error of evaluation at x
inline scalar error(const scalar x) const; inline scalar error(const scalar x) const;

View File

@ -92,6 +92,12 @@ inline Foam::scalar Foam::quadraticEqn::value(const scalar x) const
} }
inline Foam::scalar Foam::quadraticEqn::derivative(const scalar x) const
{
return x*2*a() + b();
}
inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const
{ {
return mag(SMALL*x*(x*2*a() + b())); return mag(SMALL*x*(x*2*a() + b()));

View File

@ -697,7 +697,7 @@ void Foam::motionSmootherAlgo::movePoints()
{ {
// Make sure to clear out tetPtIs since used in checks (sometimes, should // Make sure to clear out tetPtIs since used in checks (sometimes, should
// really check) // really check)
mesh_.clearAdditionalGeom(); mesh_.clearTetBasePtIs();
pp_.movePoints(mesh_.points()); pp_.movePoints(mesh_.points());
} }

View File

@ -22,11 +22,6 @@ streamLine/streamLine.C
streamLine/streamLineParticle.C streamLine/streamLineParticle.C
streamLine/streamLineParticleCloud.C streamLine/streamLineParticleCloud.C
wallBoundedStreamLine/wallBoundedStreamLine.C
wallBoundedStreamLine/wallBoundedStreamLineParticle.C
wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
wallBoundedStreamLine/wallBoundedParticle.C
surfaceInterpolate/surfaceInterpolate.C surfaceInterpolate/surfaceInterpolate.C
regionSizeDistribution/regionSizeDistribution.C regionSizeDistribution/regionSizeDistribution.C

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -30,7 +30,7 @@ License
Foam::findCellParticle::findCellParticle Foam::findCellParticle::findCellParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPtI, const label tetPtI,
@ -38,7 +38,24 @@ Foam::findCellParticle::findCellParticle
const label data const label data
) )
: :
particle(mesh, position, celli, tetFacei, tetPtI), particle(mesh, coordinates, celli, tetFacei, tetPtI),
start_(position()),
end_(end),
data_(data)
{}
Foam::findCellParticle::findCellParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const point& end,
const label data
)
:
particle(mesh, position, celli),
start_(this->position()),
end_(end), end_(end),
data_(data) data_(data)
{} {}
@ -57,15 +74,15 @@ Foam::findCellParticle::findCellParticle
{ {
if (is.format() == IOstream::ASCII) if (is.format() == IOstream::ASCII)
{ {
is >> end_; is >> start_ >> end_;
data_ = readLabel(is); data_ = readLabel(is);
} }
else else
{ {
is.read is.read
( (
reinterpret_cast<char*>(&end_), reinterpret_cast<char*>(&start_),
sizeof(end_) + sizeof(data_) sizeof(start_) + sizeof(end_) + sizeof(data_)
); );
} }
} }
@ -90,21 +107,13 @@ bool Foam::findCellParticle::move
td.switchProcessor = false; td.switchProcessor = false;
td.keepParticle = true; td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*maxTrackLen; while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
{ {
// set the lagrangian time-step const scalar f = 1 - stepFraction();
scalar dt = min(dtMax, tEnd); trackToFace(f*(end_ - start_), f, td);
dt *= trackToFace(end_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/maxTrackLen;
} }
if (tEnd < SMALL || !td.keepParticle) if (stepFraction() == 1 || !td.keepParticle)
{ {
// Hit endpoint or patch. If patch hit could do fancy stuff but just // Hit endpoint or patch. If patch hit could do fancy stuff but just
// to use the patch point is good enough for now. // to use the patch point is good enough for now.
@ -214,6 +223,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
if (os.format() == IOstream::ASCII) if (os.format() == IOstream::ASCII)
{ {
os << static_cast<const particle&>(p) os << static_cast<const particle&>(p)
<< token::SPACE << p.start_
<< token::SPACE << p.end_ << token::SPACE << p.end_
<< token::SPACE << p.data_; << token::SPACE << p.data_;
} }
@ -222,8 +232,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
os << static_cast<const particle&>(p); os << static_cast<const particle&>(p);
os.write os.write
( (
reinterpret_cast<const char*>(&p.end_), reinterpret_cast<const char*>(&p.start_),
sizeof(p.end_) + sizeof(p.data_) sizeof(p.start_) + sizeof(p.end_) + sizeof(p.data_)
); );
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -63,6 +63,9 @@ class findCellParticle
{ {
// Private data // Private data
//- Start point to track from
point start_;
//- End point to track to //- End point to track to
point end_; point end_;
@ -119,7 +122,7 @@ public:
findCellParticle findCellParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPtI, const label tetPtI,
@ -127,6 +130,17 @@ public:
const label data const label data
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology
findCellParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const point& end,
const label data
);
//- Construct from Istream //- Construct from Istream
findCellParticle findCellParticle
( (
@ -166,18 +180,42 @@ public:
// Member Functions // Member Functions
//- Point to track from
const point& start() const
{
return start_;
}
//- Point to track from
point& start()
{
return start_;
}
//- Point to track to //- Point to track to
const point& end() const const point& end() const
{ {
return end_; return end_;
} }
//- Point to track to
point& end()
{
return end_;
}
//- Transported label //- Transported label
label data() const label data() const
{ {
return data_; return data_;
} }
//- Transported label
label& data()
{
return data_;
}
// Tracking // Tracking

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -108,24 +108,16 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
const point end = start-distance_*nf[patchFacei]; const point end = start-distance_*nf[patchFacei];
// Find tet for starting location // Add a particle to the cloud with originating face as passive data
label celli = -1;
label tetFacei = -1;
label tetPtI = -1;
mesh_.findCellFacePt(start, celli, tetFacei, tetPtI);
// Add to cloud. Add originating face as passive data
cloud.addParticle cloud.addParticle
( (
new findCellParticle new findCellParticle
( (
mesh_, mesh_,
start, start,
celli, -1,
tetFacei,
tetPtI,
end, end,
globalWalls.toGlobal(nPatchFaces) // passive data globalWalls.toGlobal(nPatchFaces) // passive data
) )
); );

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -26,35 +26,8 @@ License
#include "streamLineParticle.H" #include "streamLineParticle.H"
#include "vectorFieldIOField.H" #include "vectorFieldIOField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
// defineParticleTypeNameAndDebug(streamLineParticle, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
(
trackingData& td,
const scalar dt,
const vector& U
) const
{
particle testParticle(*this);
bool oldKeepParticle = td.keepParticle;
bool oldSwitchProcessor = td.switchProcessor;
scalar fraction = testParticle.trackToFace(position()+dt*U, td);
td.keepParticle = oldKeepParticle;
td.switchProcessor = oldSwitchProcessor;
// Adapt the dt to subdivide the trajectory into substeps.
return dt*fraction/td.nSubCycle_;
}
Foam::vector Foam::streamLineParticle::interpolateFields Foam::vector Foam::streamLineParticle::interpolateFields
( (
const trackingData& td, const trackingData& td,
@ -129,7 +102,6 @@ Foam::streamLineParticle::streamLineParticle
{ {
if (readFields) if (readFields)
{ {
//if (is.format() == IOstream::ASCII)
List<scalarList> sampledScalars; List<scalarList> sampledScalars;
List<vectorList> sampledVectors; List<vectorList> sampledVectors;
@ -174,31 +146,22 @@ Foam::streamLineParticle::streamLineParticle
bool Foam::streamLineParticle::move bool Foam::streamLineParticle::move
( (
trackingData& td, trackingData& td,
const scalar trackTime const scalar
) )
{ {
streamLineParticle& p = static_cast<streamLineParticle&>(*this);
td.switchProcessor = false; td.switchProcessor = false;
td.keepParticle = true; td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*trackTime; const scalar maxDt = mesh().bounds().mag();
scalar maxDt = mesh_.bounds().mag();
while while (td.keepParticle && !td.switchProcessor && lifeTime_ > 0)
(
td.keepParticle
&& !td.switchProcessor
&& lifeTime_ > 0
)
{ {
// set the lagrangian time-step
scalar dt = maxDt; scalar dt = maxDt;
// Cross cell in steps: // Cross cell in steps:
// - at subiter 0 calculate dt to cross cell in nSubCycle steps // - at subiter 0 calculate dt to cross cell in nSubCycle steps
// - at the last subiter do all of the remaining track // - at the last subiter do all of the remaining track
for (label subIter = 0; subIter < td.nSubCycle_; subIter++) for (label subIter = 0; subIter < max(1, td.nSubCycle_); subIter++)
{ {
--lifeTime_; --lifeTime_;
@ -224,37 +187,27 @@ bool Foam::streamLineParticle::move
if (td.trackLength_ < GREAT) if (td.trackLength_ < GREAT)
{ {
// No sub-cycling. Track a set length on each step.
dt = td.trackLength_; dt = td.trackLength_;
//Pout<< " subiteration " << subIter
// << " : fixed length: updated dt:" << dt << endl;
} }
else if (subIter == 0 && td.nSubCycle_ > 1) else if (subIter == 0)
{ {
// Adapt dt to cross cell in a few steps // Sub-cycling. Cross the cell in nSubCycle steps.
dt = calcSubCycleDeltaT(td, dt, U); particle copy(*this);
copy.trackToFace(maxDt*U, 1);
dt *= (copy.stepFraction() - stepFraction())/td.nSubCycle_;
} }
else if (subIter == td.nSubCycle_ - 1) else if (subIter == td.nSubCycle_ - 1)
{ {
// Do full step on last subcycle // Sub-cycling. Track the whole cell on the last step.
dt = maxDt; dt = maxDt;
} }
trackToFace(dt*U, 0, td);
scalar fraction = trackToFace(position() + dt*U, td);
dt *= fraction;
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
if (tEnd <= ROOTVSMALL)
{
// Force removal
lifeTime_ = 0;
}
if if
( (
face() != -1 onFace()
|| !td.keepParticle || !td.keepParticle
|| td.switchProcessor || td.switchProcessor
|| lifeTime_ == 0 || lifeTime_ == 0
@ -265,17 +218,16 @@ bool Foam::streamLineParticle::move
} }
} }
if (!td.keepParticle || lifeTime_ == 0) if (!td.keepParticle || lifeTime_ == 0)
{ {
if (lifeTime_ == 0) if (lifeTime_ == 0)
{ {
// Failure exit. Particle stagnated or it's life ran out.
if (debug) if (debug)
{ {
Pout<< "streamLineParticle : Removing stagnant particle:" Pout<< "streamLineParticle: Removing stagnant particle:"
<< p.position() << position() << " sampled positions:"
<< " sampled positions:" << sampledPositions_.size() << sampledPositions_.size() << endl;
<< endl;
} }
td.keepParticle = false; td.keepParticle = false;
} }
@ -287,29 +239,25 @@ bool Foam::streamLineParticle::move
if (debug) if (debug)
{ {
Pout<< "streamLineParticle : Removing particle:" Pout<< "streamLineParticle: Removing particle:" << position()
<< p.position()
<< " sampled positions:" << sampledPositions_.size() << " sampled positions:" << sampledPositions_.size()
<< endl; << endl;
} }
} }
// Transfer particle data into trackingData. // Transfer particle data into trackingData.
//td.allPositions_.append(sampledPositions_);
td.allPositions_.append(vectorList()); td.allPositions_.append(vectorList());
vectorList& top = td.allPositions_.last(); vectorList& top = td.allPositions_.last();
top.transfer(sampledPositions_); top.transfer(sampledPositions_);
forAll(sampledScalars_, i) forAll(sampledScalars_, i)
{ {
//td.allScalars_[i].append(sampledScalars_[i]);
td.allScalars_[i].append(scalarList()); td.allScalars_[i].append(scalarList());
scalarList& top = td.allScalars_[i].last(); scalarList& top = td.allScalars_[i].last();
top.transfer(sampledScalars_[i]); top.transfer(sampledScalars_[i]);
} }
forAll(sampledVectors_, i) forAll(sampledVectors_, i)
{ {
//td.allVectors_[i].append(sampledVectors_[i]);
td.allVectors_[i].append(vectorList()); td.allVectors_[i].append(vectorList());
vectorList& top = td.allVectors_[i].last(); vectorList& top = td.allVectors_[i].last();
top.transfer(sampledVectors_[i]); top.transfer(sampledVectors_[i]);
@ -433,18 +381,11 @@ void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
); );
c.checkFieldIOobject(c, sampledPositions); c.checkFieldIOobject(c, sampledPositions);
// vectorFieldIOField sampleVelocity
// (
// c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
// );
// c.checkFieldIOobject(c, sampleVelocity);
label i = 0; label i = 0;
forAllIter(Cloud<streamLineParticle>, c, iter) forAllIter(Cloud<streamLineParticle>, c, iter)
{ {
iter().lifeTime_ = lifeTime[i]; iter().lifeTime_ = lifeTime[i];
iter().sampledPositions_.transfer(sampledPositions[i]); iter().sampledPositions_.transfer(sampledPositions[i]);
// iter().sampleVelocity_.transfer(sampleVelocity[i]);
i++; i++;
} }
} }
@ -466,24 +407,17 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
c.fieldIOobject("sampledPositions", IOobject::NO_READ), c.fieldIOobject("sampledPositions", IOobject::NO_READ),
np np
); );
// vectorFieldIOField sampleVelocity
// (
// c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
// np
// );
label i = 0; label i = 0;
forAllConstIter(Cloud<streamLineParticle>, c, iter) forAllConstIter(Cloud<streamLineParticle>, c, iter)
{ {
lifeTime[i] = iter().lifeTime_; lifeTime[i] = iter().lifeTime_;
sampledPositions[i] = iter().sampledPositions_; sampledPositions[i] = iter().sampledPositions_;
// sampleVelocity[i] = iter().sampleVelocity_;
i++; i++;
} }
lifeTime.write(); lifeTime.write();
sampledPositions.write(); sampledPositions.write();
// sampleVelocity.write();
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -46,16 +46,12 @@ SourceFiles
namespace Foam namespace Foam
{ {
class streamLineParticleCloud;
// Forward declaration of friend functions and operators // Forward declaration of friend functions and operators
class streamLineParticle; class streamLineParticle;
Ostream& operator<<(Ostream&, const streamLineParticle&); Ostream& operator<<(Ostream&, const streamLineParticle&);
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class streamLineParticle Declaration Class streamLineParticle Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -64,32 +60,38 @@ class streamLineParticle
: :
public particle public particle
{ {
public: public:
//- Class used to pass tracking data to the trackToFace function
class trackingData class trackingData
: :
public particle::TrackingData<Cloud<streamLineParticle>> public particle::TrackingData<Cloud<streamLineParticle>>
{ {
public: public:
// Public data
const PtrList<interpolation<scalar>>& vsInterp_; const PtrList<interpolation<scalar>>& vsInterp_;
const PtrList<interpolation<vector>>& vvInterp_;
const label UIndex_;
const bool trackForward_;
const label nSubCycle_;
const scalar trackLength_;
DynamicList<vectorList>& allPositions_; const PtrList<interpolation<vector>>& vvInterp_;
List<DynamicList<scalarList>>& allScalars_;
List<DynamicList<vectorList>>& allVectors_; const label UIndex_;
const bool trackForward_;
const label nSubCycle_;
const scalar trackLength_;
DynamicList<vectorList>& allPositions_;
List<DynamicList<scalarList>>& allScalars_;
List<DynamicList<vectorList>>& allVectors_;
// Constructors // Constructors
//- Construct from components
trackingData trackingData
( (
Cloud<streamLineParticle>& cloud, Cloud<streamLineParticle>& cloud,
@ -99,7 +101,6 @@ public:
const bool trackForward, const bool trackForward,
const label nSubCycle, const label nSubCycle,
const scalar trackLength, const scalar trackLength,
DynamicList<List<point>>& allPositions, DynamicList<List<point>>& allPositions,
List<DynamicList<scalarList>>& allScalars, List<DynamicList<scalarList>>& allScalars,
List<DynamicList<vectorList>>& allVectors List<DynamicList<vectorList>>& allVectors
@ -112,7 +113,6 @@ public:
trackForward_(trackForward), trackForward_(trackForward),
nSubCycle_(nSubCycle), nSubCycle_(nSubCycle),
trackLength_(trackLength), trackLength_(trackLength),
allPositions_(allPositions), allPositions_(allPositions),
allScalars_(allScalars), allScalars_(allScalars),
allVectors_(allVectors) allVectors_(allVectors)
@ -139,22 +139,6 @@ private:
// Private Member Functions // Private Member Functions
//- Estimate dt to cross from current face to next one in nSubCycle
// steps.
scalar calcSubCycleDeltaT
(
trackingData& td,
const scalar dt,
const vector& U
) const;
void constrainVelocity
(
trackingData& td,
const scalar dt,
vector& U
);
//- Interpolate all quantities; return interpolated velocity. //- Interpolate all quantities; return interpolated velocity.
vector interpolateFields vector interpolateFields
( (
@ -195,8 +179,7 @@ public:
return autoPtr<particle>(new streamLineParticle(*this)); return autoPtr<particle>(new streamLineParticle(*this));
} }
//- Factory class to read-construct particles used for //- Factory class to read-construct particles used for parallel transfer
// parallel transfer
class iNew class iNew
{ {
const polyMesh& mesh_; const polyMesh& mesh_;
@ -223,8 +206,7 @@ public:
// Tracking // Tracking
//- Track all particles to their end point //- Track all particles to their end point
bool move(trackingData&, const scalar trackTime); bool move(trackingData&, const scalar);
//- Overridable function to handle the particle hitting a patch //- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions // Executed before other patch-hitting functions

View File

@ -1,425 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "wallBoundedParticle.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const std::size_t Foam::wallBoundedParticle::sizeofFields_
(
sizeof(wallBoundedParticle) - sizeof(particle)
);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::edge Foam::wallBoundedParticle::currentEdge() const
{
if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
{
FatalErrorInFunction
<< "Particle:"
<< info()
<< "cannot both be on a mesh edge and a face-diagonal edge."
<< " meshEdgeStart_:" << meshEdgeStart_
<< " diagEdge_:" << diagEdge_
<< abort(FatalError);
}
const Foam::face& f = mesh_.faces()[tetFace()];
if (meshEdgeStart_ != -1)
{
return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
}
else
{
label faceBasePtI = mesh_.tetBasePtIs()[tetFace()];
label diagPtI = (faceBasePtI+diagEdge_)%f.size();
return edge(f[faceBasePtI], f[diagPtI]);
}
}
void Foam::wallBoundedParticle::crossEdgeConnectedFace
(
const edge& meshEdge
)
{
// Update tetFace, tetPt
particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
// Update face to be same as tracking one
face() = tetFace();
// And adapt meshEdgeStart_.
const Foam::face& f = mesh_.faces()[tetFace()];
label fp = findIndex(f, meshEdge[0]);
if (f.nextLabel(fp) == meshEdge[1])
{
meshEdgeStart_ = fp;
}
else
{
label fpMin1 = f.rcIndex(fp);
if (f[fpMin1] == meshEdge[1])
{
meshEdgeStart_ = fpMin1;
}
else
{
FatalErrorInFunction
<< "Problem :"
<< " particle:"
<< info()
<< "face:" << tetFace()
<< " verts:" << f
<< " meshEdge:" << meshEdge
<< abort(FatalError);
}
}
diagEdge_ = -1;
// Check that still on same mesh edge
const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
if (eNew != meshEdge)
{
FatalErrorInFunction
<< "Problem" << abort(FatalError);
}
}
void Foam::wallBoundedParticle::crossDiagonalEdge()
{
if (diagEdge_ == -1)
{
FatalErrorInFunction
<< "Particle:"
<< info()
<< "not on a diagonal edge" << abort(FatalError);
}
if (meshEdgeStart_ != -1)
{
FatalErrorInFunction
<< "Particle:"
<< info()
<< "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
}
const Foam::face& f = mesh_.faces()[tetFace()];
// tetPtI starts from 1, goes up to f.size()-2
if (tetPt() == diagEdge_)
{
tetPt() = f.rcIndex(tetPt());
}
else
{
label nextTetPt = f.fcIndex(tetPt());
if (diagEdge_ == nextTetPt)
{
tetPt() = nextTetPt;
}
else
{
FatalErrorInFunction
<< "Particle:"
<< info()
<< "tetPt:" << tetPt()
<< " diagEdge:" << diagEdge_ << abort(FatalError);
}
}
meshEdgeStart_ = -1;
}
Foam::scalar Foam::wallBoundedParticle::trackFaceTri
(
const vector& endPosition,
label& minEdgeI
)
{
// Track p from position to endPosition
const triFace tri(currentTetIndices().faceTriIs(mesh_));
vector n = tri.normal(mesh_.points());
n /= mag(n)+VSMALL;
// Check which edge intersects the trajectory.
// Project trajectory onto triangle
minEdgeI = -1;
scalar minS = 1; // end position
edge currentE(-1, -1);
if (meshEdgeStart_ != -1 || diagEdge_ != -1)
{
currentE = currentEdge();
}
// Determine path along line position+s*d to see where intersections
// are.
forAll(tri, i)
{
label j = tri.fcIndex(i);
const point& pt0 = mesh_.points()[tri[i]];
const point& pt1 = mesh_.points()[tri[j]];
if (edge(tri[i], tri[j]) == currentE)
{
// Do not check particle is on
continue;
}
// Outwards pointing normal
vector edgeNormal = (pt1-pt0)^n;
edgeNormal /= mag(edgeNormal)+VSMALL;
// Determine whether position and end point on either side of edge.
scalar sEnd = (endPosition-pt0)&edgeNormal;
if (sEnd >= 0)
{
// endPos is outside triangle. position() should always be
// inside.
scalar sStart = (position()-pt0)&edgeNormal;
if (mag(sEnd-sStart) > VSMALL)
{
scalar s = sStart/(sStart-sEnd);
if (s >= 0 && s < minS)
{
minS = s;
minEdgeI = i;
}
}
}
}
if (minEdgeI != -1)
{
position() += minS*(endPosition-position());
}
else
{
// Did not hit any edge so tracked to the end position. Set position
// without any calculation to avoid truncation errors.
position() = endPosition;
minS = 1.0;
}
// Project position() back onto plane of triangle
const point& triPt = mesh_.points()[tri[0]];
position() -= ((position()-triPt)&n)*n;
return minS;
}
bool Foam::wallBoundedParticle::isTriAlongTrack
(
const point& endPosition
) const
{
const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
const edge currentE = currentEdge();
//if (debug)
{
if
(
currentE[0] == currentE[1]
|| findIndex(triVerts, currentE[0]) == -1
|| findIndex(triVerts, currentE[1]) == -1
)
{
FatalErrorInFunction
<< "Edge " << currentE << " not on triangle " << triVerts
<< info()
<< abort(FatalError);
}
}
const vector dir = endPosition-position();
// Get normal of currentE
vector n = triVerts.normal(mesh_.points());
n /= mag(n);
forAll(triVerts, i)
{
label j = triVerts.fcIndex(i);
const point& pt0 = mesh_.points()[triVerts[i]];
const point& pt1 = mesh_.points()[triVerts[j]];
if (edge(triVerts[i], triVerts[j]) == currentE)
{
vector edgeNormal = (pt1-pt0)^n;
return (dir&edgeNormal) < 0;
}
}
FatalErrorInFunction
<< "Problem" << abort(FatalError);
return false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::wallBoundedParticle::wallBoundedParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const label tetFacei,
const label tetPtI,
const label meshEdgeStart,
const label diagEdge
)
:
particle(mesh, position, celli, tetFacei, tetPtI),
meshEdgeStart_(meshEdgeStart),
diagEdge_(diagEdge)
{}
Foam::wallBoundedParticle::wallBoundedParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
)
:
particle(mesh, is, readFields)
{
if (readFields)
{
if (is.format() == IOstream::ASCII)
{
is >> meshEdgeStart_ >> diagEdge_;
}
else
{
is.read
(
reinterpret_cast<char*>(&meshEdgeStart_),
sizeof(meshEdgeStart_)
+ sizeof(diagEdge_)
);
}
}
// Check state of Istream
is.check
(
"wallBoundedParticle::wallBoundedParticle"
"(const Cloud<wallBoundedParticle>&, Istream&, bool)"
);
}
Foam::wallBoundedParticle::wallBoundedParticle
(
const wallBoundedParticle& p
)
:
particle(p),
meshEdgeStart_(p.meshEdgeStart_),
diagEdge_(p.diagEdge_)
{}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const wallBoundedParticle& p
)
{
if (os.format() == IOstream::ASCII)
{
os << static_cast<const particle&>(p)
<< token::SPACE << p.meshEdgeStart_
<< token::SPACE << p.diagEdge_;
}
else
{
os << static_cast<const particle&>(p);
os.write
(
reinterpret_cast<const char*>(&p.meshEdgeStart_),
wallBoundedParticle::sizeofFields_
);
}
return os;
}
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const InfoProxy<wallBoundedParticle>& ip
)
{
const wallBoundedParticle& p = ip.t_;
tetPointRef tpr(p.currentTet());
os << " " << static_cast<const particle&>(p) << nl
<< " tet:" << nl;
os << " ";
meshTools::writeOBJ(os, tpr.a());
os << " ";
meshTools::writeOBJ(os, tpr.b());
os << " ";
meshTools::writeOBJ(os, tpr.c());
os << " ";
meshTools::writeOBJ(os, tpr.d());
os << " l 1 2" << nl
<< " l 1 3" << nl
<< " l 1 4" << nl
<< " l 2 3" << nl
<< " l 2 4" << nl
<< " l 3 4" << nl;
os << " ";
meshTools::writeOBJ(os, p.position());
return os;
}
// ************************************************************************* //

View File

@ -1,364 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
Class
Foam::wallBoundedParticle
Description
Particle class that tracks on triangles of boundary faces. Use
trackToEdge similar to trackToFace on particle.
SourceFiles
wallBoundedParticle.C
wallBoundedParticleTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef wallBoundedParticle_H
#define wallBoundedParticle_H
#include "particle.H"
#include "autoPtr.H"
#include "InfoProxy.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
class wallBoundedParticle;
Ostream& operator<<(Ostream&, const wallBoundedParticle&);
Ostream& operator<<(Ostream&, const InfoProxy<wallBoundedParticle>&);
/*---------------------------------------------------------------------------*\
Class wallBoundedParticle Declaration
\*---------------------------------------------------------------------------*/
class wallBoundedParticle
:
public particle
{
// Private data
//- Size in bytes of the fields
static const std::size_t sizeofFields_;
public:
//- Class used to pass tracking data to the trackToFace function
template<class CloudType>
class TrackingData
:
public particle::TrackingData<CloudType>
{
public:
const PackedBoolList& isWallPatch_;
// Constructors
inline TrackingData
(
CloudType& cloud,
const PackedBoolList& isWallPatch
)
:
particle::TrackingData<CloudType>
(
cloud
),
isWallPatch_(isWallPatch)
{}
};
protected:
// Protected data
//- Particle is on mesh edge:
// const face& f = mesh.faces()[tetFace()]
// const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
// Note that this real edge
// is also one of the edges of the face-triangle (from
// tetFace()+tetPt()).
label meshEdgeStart_;
//- Particle is on diagonal edge:
// const face& f = mesh.faces()[tetFace()]
// label faceBasePtI = mesh.tetBasePtIs()[facei];
// label diagPtI = (faceBasePtI+diagEdge_)%f.size();
// const edge e(f[faceBasePtI], f[diagPtI]);
label diagEdge_;
// Protected Member Functions
//- Construct current edge
edge currentEdge() const;
//- Check if inside current tet
//void checkInside() const;
//- Check if on current edge
//void checkOnEdge() const;
//- Check if point on triangle
//void checkOnTriangle(const point&) const;
//- Cross mesh edge into different face on same cell
void crossEdgeConnectedFace(const edge& meshEdge);
//- Cross diagonal edge into different triangle on same face,cell
void crossDiagonalEdge();
//- Track through single triangle
scalar trackFaceTri(const vector& endPosition, label& minEdgeI);
//- Is current triangle in the track direction
bool isTriAlongTrack(const point& endPosition) const;
// Patch interactions
//- Do all patch interaction
template<class TrackData>
void patchInteraction(TrackData& td, const scalar trackFraction);
//- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions
template<class TrackData>
bool hitPatch
(
const polyPatch&,
TrackData& td,
const label patchi,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a wedge
template<class TrackData>
void hitWedgePatch
(
const wedgePolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
// symmetry plane
template<class TrackData>
void hitSymmetryPlanePatch
(
const symmetryPlanePolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
// symmetry patch
template<class TrackData>
void hitSymmetryPatch
(
const symmetryPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a cyclic
template<class TrackData>
void hitCyclicPatch
(
const cyclicPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
//- processorPatch
template<class TrackData>
void hitProcessorPatch
(
const processorPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a wallPatch
template<class TrackData>
void hitWallPatch
(
const wallPolyPatch&,
TrackData& td,
const tetIndices&
);
//- Overridable function to handle the particle hitting a polyPatch
template<class TrackData>
void hitPatch
(
const polyPatch&,
TrackData& td
);
public:
// Constructors
//- Construct from components
wallBoundedParticle
(
const polyMesh& c,
const vector& position,
const label celli,
const label tetFacei,
const label tetPtI,
const label meshEdgeStart,
const label diagEdge
);
//- Construct from Istream
wallBoundedParticle
(
const polyMesh& c,
Istream& is,
bool readFields = true
);
//- Construct copy
wallBoundedParticle(const wallBoundedParticle& p);
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new wallBoundedParticle(*this));
}
//- Factory class to read-construct particles used for
// parallel transfer
class iNew
{
const polyMesh& mesh_;
public:
iNew(const polyMesh& mesh)
:
mesh_(mesh)
{}
autoPtr<wallBoundedParticle> operator()
(
Istream& is
) const
{
return autoPtr<wallBoundedParticle>
(
new wallBoundedParticle(mesh_, is, true)
);
}
};
// Member Functions
// Access
//- -1 or label of mesh edge
inline label meshEdgeStart() const
{
return meshEdgeStart_;
}
//- -1 or diagonal edge
inline label diagEdge() const
{
return diagEdge_;
}
// Track
//- Equivalent of trackToFace
template<class TrackData>
scalar trackToEdge
(
TrackData& td,
const vector& endPosition
);
// Info
//- Return info proxy.
// Used to print particle information to a stream
inline InfoProxy<wallBoundedParticle> info() const
{
return *this;
}
// I-O
//- Read
template<class CloudType>
static void readFields(CloudType&);
//- Write
template<class CloudType>
static void writeFields(const CloudType&);
// Ostream Operator
friend Ostream& operator<<
(
Ostream&,
const wallBoundedParticle&
);
friend Ostream& operator<<
(
Ostream&,
const InfoProxy<wallBoundedParticle>&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "wallBoundedParticleTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,559 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "wallBoundedParticle.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class TrackData>
void Foam::wallBoundedParticle::patchInteraction
(
TrackData& td,
const scalar trackFraction
)
{
typedef typename TrackData::cloudType::particleType particleType;
particleType& p = static_cast<particleType&>(*this);
p.hitFace(td);
if (!internalFace(facei_))
{
label origFacei = facei_;
label patchi = patch(facei_);
// No action taken for tetPti_ for tetFacei_ here, handled by
// patch interaction call or later during processor transfer.
// Dummy tet indices. What to do here?
tetIndices faceHitTetIs;
if
(
!p.hitPatch
(
mesh_.boundaryMesh()[patchi],
td,
patchi,
trackFraction,
faceHitTetIs
)
)
{
// Did patch interaction model switch patches?
// Note: recalculate meshEdgeStart_, diagEdge_!
if (facei_ != origFacei)
{
patchi = patch(facei_);
}
const polyPatch& patch = mesh_.boundaryMesh()[patchi];
if (isA<wedgePolyPatch>(patch))
{
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPlanePolyPatch>(patch))
{
p.hitSymmetryPlanePatch
(
static_cast<const symmetryPlanePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
);
}
else
{
p.hitPatch(patch, td);
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class TrackData>
Foam::scalar Foam::wallBoundedParticle::trackToEdge
(
TrackData& td,
const vector& endPosition
)
{
// Track particle to a given position and returns 1.0 if the
// trajectory is completed without hitting a face otherwise
// stops at the face and returns the fraction of the trajectory
// completed.
// on entry 'stepFraction()' should be set to the fraction of the
// time-step at which the tracking starts.
// Are we on a track face? If not we do a topological walk.
// Particle:
// - cell_ always set
// - tetFace_, tetPt_ always set (these identify tet particle is in)
// - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
//checkInside();
//checkOnTriangle(position());
//if (meshEdgeStart_ != -1 || diagEdge_ != -1)
//{
// checkOnEdge();
//}
scalar trackFraction = 0.0;
if (!td.isWallPatch_[tetFace()])
{
// Don't track across face. Just walk in cell. Particle is on
// mesh edge (as indicated by meshEdgeStart_).
const edge meshEdge(currentEdge());
// If internal face check whether to go to neighbour cell or just
// check to the other internal tet on the edge.
if (mesh_.isInternalFace(tetFace()))
{
label nbrCelli =
(
celli_ == mesh_.faceOwner()[facei_]
? mesh_.faceNeighbour()[facei_]
: mesh_.faceOwner()[facei_]
);
// Check angle to nbrCell tet. Is it in the direction of the
// endposition? I.e. since volume of nbr tet is positive the
// tracking direction should be into the tet.
tetIndices nbrTi(nbrCelli, tetFacei_, tetPti_, mesh_);
if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
{
// Change into nbrCell. No need to change tetFace, tetPt.
//Pout<< " crossed from cell:" << celli_
// << " into " << nbrCelli << endl;
celli_ = nbrCelli;
patchInteraction(td, trackFraction);
}
else
{
// Walk to other face on edge. Changes tetFace, tetPt but not
// cell.
crossEdgeConnectedFace(meshEdge);
patchInteraction(td, trackFraction);
}
}
else
{
// Walk to other face on edge. This might give loop since
// particle should have been removed?
crossEdgeConnectedFace(meshEdge);
patchInteraction(td, trackFraction);
}
}
else
{
// We're inside a tet on the wall. Check if the current tet is
// the one to cross. If not we cross into the neighbouring triangle.
if (mesh_.isInternalFace(tetFace()))
{
FatalErrorInFunction
<< "Can only track on boundary faces."
<< " Face:" << tetFace()
<< " at:" << mesh_.faceCentres()[tetFace()]
<< abort(FatalError);
}
point projectedEndPosition = endPosition;
// Remove normal component
{
const triFace tri(currentTetIndices().faceTriIs(mesh_));
vector n = tri.normal(mesh_.points());
n /= mag(n);
const point& basePt = mesh_.points()[tri[0]];
projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
}
bool doTrack = false;
if (meshEdgeStart_ == -1 && diagEdge_ == -1)
{
// We're starting and not yet on an edge.
doTrack = true;
}
else
{
// See if the current triangle has got a point on the
// correct side of the edge.
doTrack = isTriAlongTrack(projectedEndPosition);
}
if (doTrack)
{
// Track across triangle. Return triangle edge crossed.
label triEdgeI = -1;
trackFraction = trackFaceTri(projectedEndPosition, triEdgeI);
if (triEdgeI == -1)
{
// Reached endpoint
//checkInside();
diagEdge_ = -1;
meshEdgeStart_ = -1;
return trackFraction;
}
const tetIndices ti(currentTetIndices());
// Triangle (faceTriIs) gets constructed from
// f[faceBasePtI_],
// f[facePtAI_],
// f[facePtBI_]
//
// So edge indices are:
// 0 : edge between faceBasePtI_ and facePtAI_
// 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
// 2 : edge between facePtBI_ and faceBasePtI_
const Foam::face& f = mesh_.faces()[ti.face()];
const label fp0 = ti.faceBasePt();
if (triEdgeI == 0)
{
if (ti.facePtA() == f.fcIndex(fp0))
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = fp0;
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
}
else if (ti.facePtA() == f.rcIndex(fp0))
{
//Note: should not happen since boundary face so owner
//Pout<< "Real edge." << endl;
FatalErrorInFunction
<< abort(FatalError);
diagEdge_ = -1;
meshEdgeStart_ = f.rcIndex(fp0);
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
}
else
{
// Get index of triangle on other side of edge.
diagEdge_ = ti.facePtA()-fp0;
if (diagEdge_ < 0)
{
diagEdge_ += f.size();
}
meshEdgeStart_ = -1;
//checkOnEdge();
crossDiagonalEdge();
}
}
else if (triEdgeI == 1)
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = ti.facePtA();
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
}
else // if (triEdgeI == 2)
{
if (ti.facePtB() == f.rcIndex(fp0))
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = ti.facePtB();
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
}
else if (ti.facePtB() == f.fcIndex(fp0))
{
//Note: should not happen since boundary face so owner
//Pout<< "Real edge." << endl;
FatalErrorInFunction
<< abort(FatalError);
diagEdge_ = -1;
meshEdgeStart_ = fp0;
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
}
else
{
//Pout<< "Triangle edge." << endl;
// Get index of triangle on other side of edge.
diagEdge_ = ti.facePtB()-fp0;
if (diagEdge_ < 0)
{
diagEdge_ += f.size();
}
meshEdgeStart_ = -1;
//checkOnEdge();
crossDiagonalEdge();
}
}
}
else
{
// Current tet is not the right one. Check the neighbour tet.
if (meshEdgeStart_ != -1)
{
// Particle is on mesh edge so change into other face on cell
crossEdgeConnectedFace(currentEdge());
//checkOnEdge();
patchInteraction(td, trackFraction);
}
else
{
// Particle is on diagonal edge so change into the other
// triangle.
crossDiagonalEdge();
//checkOnEdge();
}
}
}
//checkInside();
return trackFraction;
}
template<class TrackData>
bool Foam::wallBoundedParticle::hitPatch
(
const polyPatch&,
TrackData& td,
const label patchi,
const scalar trackFraction,
const tetIndices& tetIs
)
{
// Disable generic patch interaction
return false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitWedgePatch
(
const wedgePolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitSymmetryPlanePatch
(
const symmetryPlanePolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitSymmetryPatch
(
const symmetryPolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitCyclicPatch
(
const cyclicPolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitProcessorPatch
(
const processorPolyPatch& pp,
TrackData& td
)
{
// Switch particle
td.switchProcessor = true;
// Adapt edgeStart_ for other side.
// E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
// on the other side between 2 and 3 so edgeStart_ should be
// f.size()-edgeStart_-1.
const Foam::face& f = mesh_.faces()[face()];
if (meshEdgeStart_ != -1)
{
meshEdgeStart_ = f.size()-meshEdgeStart_-1;
}
else
{
// diagEdge_ is relative to faceBasePt
diagEdge_ = f.size()-diagEdge_;
}
}
template<class TrackData>
void Foam::wallBoundedParticle::hitWallPatch
(
const wallPolyPatch& wpp,
TrackData& td,
const tetIndices&
)
{}
template<class TrackData>
void Foam::wallBoundedParticle::hitPatch
(
const polyPatch& wpp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class CloudType>
void Foam::wallBoundedParticle::readFields(CloudType& c)
{
if (!c.size())
{
return;
}
particle::readFields(c);
IOField<label> meshEdgeStart
(
c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
);
IOField<label> diagEdge
(
c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, diagEdge);
label i = 0;
forAllIter(typename CloudType, c, iter)
{
iter().meshEdgeStart_ = meshEdgeStart[i];
iter().diagEdge_ = diagEdge[i];
i++;
}
}
template<class CloudType>
void Foam::wallBoundedParticle::writeFields(const CloudType& c)
{
particle::writeFields(c);
label np = c.size();
IOField<label> meshEdgeStart
(
c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
np
);
IOField<label> diagEdge
(
c.fieldIOobject("diagEdge", IOobject::NO_READ),
np
);
label i = 0;
forAllConstIter(typename CloudType, c, iter)
{
meshEdgeStart[i] = iter().meshEdgeStart_;
diagEdge[i] = iter().diagEdge_;
i++;
}
meshEdgeStart.write();
diagEdge.write();
}
// ************************************************************************* //

View File

@ -1,824 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Pstream.H"
#include "functionObjectList.H"
#include "wallBoundedStreamLine.H"
#include "fvMesh.H"
#include "wallBoundedStreamLineParticleCloud.H"
#include "ReadFields.H"
#include "meshSearch.H"
#include "sampledSet.H"
#include "globalIndex.H"
#include "mapDistribute.H"
#include "interpolationCellPoint.H"
#include "PatchTools.H"
#include "meshSearchMeshObject.H"
#include "faceSet.H"
#include "mapPolyMesh.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
defineTypeNameAndDebug(wallBoundedStreamLine, 0);
addToRunTimeSelectionTable
(
functionObject,
wallBoundedStreamLine,
dictionary
);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::functionObjects::wallBoundedStreamLine::wallPatch() const
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
label nFaces = 0;
forAll(patches, patchi)
{
//if (!polyPatch::constraintType(patches[patchi].type()))
if (isA<wallPolyPatch>(patches[patchi]))
{
nFaces += patches[patchi].size();
}
}
labelList addressing(nFaces);
nFaces = 0;
forAll(patches, patchi)
{
//if (!polyPatch::constraintType(patches[patchi].type()))
if (isA<wallPolyPatch>(patches[patchi]))
{
const polyPatch& pp = patches[patchi];
forAll(pp, i)
{
addressing[nFaces++] = pp.start()+i;
}
}
}
return autoPtr<indirectPrimitivePatch>
(
new indirectPrimitivePatch
(
IndirectList<face>
(
mesh_.faces(),
addressing
),
mesh_.points()
)
);
}
Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
(
const PackedBoolList& isWallPatch,
const point& seedPt,
const label celli
) const
{
const cell& cFaces = mesh_.cells()[celli];
label minFacei = -1;
label minTetPtI = -1;
scalar minDistSqr = sqr(GREAT);
forAll(cFaces, cFacei)
{
label facei = cFaces[cFacei];
if (isWallPatch[facei])
{
const face& f = mesh_.faces()[facei];
const label fp0 = mesh_.tetBasePtIs()[facei];
const point& basePoint = mesh_.points()[f[fp0]];
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); i++)
{
const point& thisPoint = mesh_.points()[f[fp]];
label nextFp = f.fcIndex(fp);
const point& nextPoint = mesh_.points()[f[nextFp]];
const triPointRef tri(basePoint, thisPoint, nextPoint);
scalar d2 = magSqr(tri.centre() - seedPt);
if (d2 < minDistSqr)
{
minDistSqr = d2;
minFacei = facei;
minTetPtI = i-1;
}
fp = nextFp;
}
}
}
// Put particle in tet
return tetIndices
(
celli,
minFacei,
minTetPtI,
mesh_
);
}
void Foam::functionObjects::wallBoundedStreamLine::track()
{
// Determine the 'wall' patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// These are the faces that need to be followed
autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
PackedBoolList isWallPatch(mesh_.nFaces());
forAll(boundaryPatch().addressing(), i)
{
isWallPatch[boundaryPatch().addressing()[i]] = 1;
}
// Find nearest wall particle for the seedPoints
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
IDLList<wallBoundedStreamLineParticle> initialParticles;
wallBoundedStreamLineParticleCloud particles
(
mesh_,
cloudName_,
initialParticles
);
{
// Get the seed points
// ~~~~~~~~~~~~~~~~~~~
const sampledSet& seedPoints = sampledSetPtr_();
forAll(seedPoints, i)
{
const point& seedPt = seedPoints[i];
label celli = seedPoints.cells()[i];
tetIndices ids(findNearestTet(isWallPatch, seedPt, celli));
if (ids.face() != -1 && isWallPatch[ids.face()])
{
//Pout<< "Seeding particle :" << nl
// << " seedPt:" << seedPt << nl
// << " face :" << ids.face() << nl
// << " at :" << mesh_.faceCentres()[ids.face()] << nl
// << " cell :" << mesh_.cellCentres()[ids.cell()] << nl
// << endl;
particles.addParticle
(
new wallBoundedStreamLineParticle
(
mesh_,
ids.faceTri(mesh_).centre(),
ids.cell(),
ids.face(), // tetFace
ids.tetPt(),
-1, // not on a mesh edge
-1, // not on a diagonal edge
lifeTime_ // lifetime
)
);
}
else
{
Pout<< type() << " : ignoring seed " << seedPt
<< " since not in wall cell." << endl;
}
}
}
label nSeeds = returnReduce(particles.size(), sumOp<label>());
Info<< type() << " : seeded " << nSeeds << " particles." << endl;
// Read or lookup fields
PtrList<volScalarField> vsFlds;
PtrList<interpolation<scalar>> vsInterp;
PtrList<volVectorField> vvFlds;
PtrList<interpolation<vector>> vvInterp;
label UIndex = -1;
label nScalar = 0;
label nVector = 0;
forAll(fields_, i)
{
if (mesh_.foundObject<volScalarField>(fields_[i]))
{
nScalar++;
}
else if (mesh_.foundObject<volVectorField>(fields_[i]))
{
nVector++;
}
else
{
FatalErrorInFunction
<< "Cannot find field " << fields_[i] << endl
<< "Valid scalar fields are:"
<< mesh_.names(volScalarField::typeName) << endl
<< "Valid vector fields are:"
<< mesh_.names(volVectorField::typeName)
<< exit(FatalError);
}
}
vsInterp.setSize(nScalar);
nScalar = 0;
vvInterp.setSize(nVector);
nVector = 0;
forAll(fields_, i)
{
if (mesh_.foundObject<volScalarField>(fields_[i]))
{
const volScalarField& f = mesh_.lookupObject<volScalarField>
(
fields_[i]
);
vsInterp.set
(
nScalar++,
interpolation<scalar>::New
(
interpolationScheme_,
f
)
);
}
else if (mesh_.foundObject<volVectorField>(fields_[i]))
{
const volVectorField& f = mesh_.lookupObject<volVectorField>
(
fields_[i]
);
if (f.name() == UName_)
{
UIndex = nVector;
}
vvInterp.set
(
nVector++,
interpolation<vector>::New
(
interpolationScheme_,
f
)
);
}
}
// Store the names
scalarNames_.setSize(vsInterp.size());
forAll(vsInterp, i)
{
scalarNames_[i] = vsInterp[i].psi().name();
}
vectorNames_.setSize(vvInterp.size());
forAll(vvInterp, i)
{
vectorNames_[i] = vvInterp[i].psi().name();
}
// Check that we know the index of U in the interpolators.
if (UIndex == -1)
{
FatalErrorInFunction
<< "Cannot find field to move particles with : " << UName_
<< endl
<< "This field has to be present in the sampled fields "
<< fields_
<< " and in the objectRegistry." << endl
<< exit(FatalError);
}
// Sampled data
// ~~~~~~~~~~~~
// Size to maximum expected sizes.
allTracks_.clear();
allTracks_.setCapacity(nSeeds);
allScalars_.setSize(vsInterp.size());
forAll(allScalars_, i)
{
allScalars_[i].clear();
allScalars_[i].setCapacity(nSeeds);
}
allVectors_.setSize(vvInterp.size());
forAll(allVectors_, i)
{
allVectors_[i].clear();
allVectors_[i].setCapacity(nSeeds);
}
// additional particle info
wallBoundedStreamLineParticle::trackingData td
(
particles,
vsInterp,
vvInterp,
UIndex, // index of U in vvInterp
trackForward_, // track in +u direction?
trackLength_, // fixed track length
isWallPatch, // which faces are to follow
allTracks_,
allScalars_,
allVectors_
);
// Set very large dt. Note: cannot use GREAT since 1/GREAT is SMALL
// which is a trigger value for the tracking...
const scalar trackTime = Foam::sqrt(GREAT);
// Track
particles.move(td, trackTime);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::wallBoundedStreamLine::wallBoundedStreamLine
(
const word& name,
const Time& runTime,
const dictionary& dict
)
:
fvMeshFunctionObject(name, runTime, dict),
dict_(dict)
{
read(dict_);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::functionObjects::wallBoundedStreamLine::~wallBoundedStreamLine()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::functionObjects::wallBoundedStreamLine::read(const dictionary& dict)
{
if (dict != dict_)
{
dict_ = dict;
}
Info<< type() << " " << name() << ":" << nl;
dict.lookup("fields") >> fields_;
dict.lookup("U") >> UName_;
if (findIndex(fields_, UName_) == -1)
{
FatalIOErrorInFunction
(
dict
) << "Velocity field for tracking " << UName_
<< " should be present in the list of fields " << fields_
<< exit(FatalIOError);
}
dict.lookup("trackForward") >> trackForward_;
dict.lookup("lifeTime") >> lifeTime_;
if (lifeTime_ < 1)
{
FatalErrorInFunction
<< "Illegal value " << lifeTime_ << " for lifeTime"
<< exit(FatalError);
}
trackLength_ = VGREAT;
if (dict.found("trackLength"))
{
dict.lookup("trackLength") >> trackLength_;
Info<< type() << " : fixed track length specified : "
<< trackLength_ << nl << endl;
}
interpolationScheme_ = dict.lookupOrDefault
(
"interpolationScheme",
interpolationCellPoint<scalar>::typeName
);
cloudName_ = dict.lookupOrDefault<word>
(
"cloudName",
"wallBoundedStreamLine"
);
sampledSetPtr_ = sampledSet::New
(
"seedSampleSet",
mesh_,
meshSearchMeshObject::New(mesh_),
dict.subDict("seedSampleSet")
);
sampledSetAxis_ = sampledSetPtr_->axis();
scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
// Make sure that the mesh is trackable
if (debug)
{
// 1. positive volume decomposition tets
faceSet faces(mesh_, "lowQualityTetFaces", mesh_.nFaces()/100+1);
if
(
polyMeshTetDecomposition::checkFaceTets
(
mesh_,
polyMeshTetDecomposition::minTetQuality,
true,
&faces
)
)
{
label nFaces = returnReduce(faces.size(), sumOp<label>());
WarningInFunction
<< "Found " << nFaces
<<" faces with low quality or negative volume "
<< "decomposition tets. Writing to faceSet " << faces.name()
<< endl;
}
// 2. all edges on a cell having two faces
EdgeMap<label> numFacesPerEdge;
forAll(mesh_.cells(), celli)
{
const cell& cFaces = mesh_.cells()[celli];
numFacesPerEdge.clear();
forAll(cFaces, cFacei)
{
label facei = cFaces[cFacei];
const face& f = mesh_.faces()[facei];
forAll(f, fp)
{
const edge e(f[fp], f.nextLabel(fp));
EdgeMap<label>::iterator eFnd =
numFacesPerEdge.find(e);
if (eFnd != numFacesPerEdge.end())
{
eFnd()++;
}
else
{
numFacesPerEdge.insert(e, 1);
}
}
}
forAllConstIter(EdgeMap<label>, numFacesPerEdge, iter)
{
if (iter() != 2)
{
FatalErrorInFunction
<< "problem cell:" << celli
<< abort(FatalError);
}
}
}
}
return true;
}
bool Foam::functionObjects::wallBoundedStreamLine::execute()
{
return true;
}
bool Foam::functionObjects::wallBoundedStreamLine::write()
{
const Time& runTime = obr_.time();
// Do all injection and tracking
track();
if (Pstream::parRun())
{
// Append slave tracks to master ones
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
globalIndex globalTrackIDs(allTracks_.size());
// Construct a distribution map to pull all to the master.
labelListList sendMap(Pstream::nProcs());
labelListList recvMap(Pstream::nProcs());
if (Pstream::master())
{
// Master: receive all. My own first, then consecutive
// processors.
label trackI = 0;
forAll(recvMap, proci)
{
labelList& fromProc = recvMap[proci];
fromProc.setSize(globalTrackIDs.localSize(proci));
forAll(fromProc, i)
{
fromProc[i] = trackI++;
}
}
}
labelList& toMaster = sendMap[0];
toMaster.setSize(globalTrackIDs.localSize());
forAll(toMaster, i)
{
toMaster[i] = i;
}
const mapDistribute distMap
(
globalTrackIDs.size(),
sendMap.xfer(),
recvMap.xfer()
);
// Distribute the track positions. Note: use scheduled comms
// to prevent buffering.
allTracks_.shrink();
mapDistributeBase::distribute
(
Pstream::commsTypes::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allTracks_,
flipOp()
);
// Distribute the scalars
forAll(allScalars_, scalarI)
{
allScalars_[scalarI].shrink();
mapDistributeBase::distribute
(
Pstream::commsTypes::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allScalars_[scalarI],
flipOp()
);
allScalars_[scalarI].setCapacity(allScalars_[scalarI].size());
}
// Distribute the vectors
forAll(allVectors_, vectorI)
{
allVectors_[vectorI].shrink();
mapDistributeBase::distribute
(
Pstream::commsTypes::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allVectors_[vectorI],
flipOp()
);
allVectors_[vectorI].setCapacity(allVectors_[vectorI].size());
}
}
label n = 0;
forAll(allTracks_, trackI)
{
n += allTracks_[trackI].size();
}
Info<< " Tracks:" << allTracks_.size() << nl
<< " Total samples:" << n << endl;
// Massage into form suitable for writers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (Pstream::master() && allTracks_.size())
{
// Make output directory
fileName vtkPath
(
Pstream::parRun()
? runTime.path()/".."/"postProcessing"/"sets"/name()
: runTime.path()/"postProcessing"/"sets"/name()
);
if (mesh_.name() != fvMesh::defaultRegion)
{
vtkPath = vtkPath/mesh_.name();
}
vtkPath = vtkPath/mesh_.time().timeName();
mkDir(vtkPath);
// Convert track positions
PtrList<coordSet> tracks(allTracks_.size());
forAll(allTracks_, trackI)
{
tracks.set
(
trackI,
new coordSet
(
"track" + Foam::name(trackI),
sampledSetAxis_ //"xyz"
)
);
tracks[trackI].transfer(allTracks_[trackI]);
}
// Convert scalar values
if (allScalars_.size() > 0)
{
List<List<scalarField>> scalarValues(allScalars_.size());
forAll(allScalars_, scalarI)
{
DynamicList<scalarList>& allTrackVals =
allScalars_[scalarI];
scalarValues[scalarI].setSize(allTrackVals.size());
forAll(allTrackVals, trackI)
{
scalarList& trackVals = allTrackVals[trackI];
scalarValues[scalarI][trackI].transfer(trackVals);
}
}
fileName vtkFile
(
vtkPath
/ scalarFormatterPtr_().getFileName
(
tracks[0],
scalarNames_
)
);
Info<< " Writing data to " << vtkFile.path() << endl;
scalarFormatterPtr_().write
(
true, // writeTracks
tracks,
scalarNames_,
scalarValues,
OFstream(vtkFile)()
);
}
// Convert vector values
if (allVectors_.size() > 0)
{
List<List<vectorField>> vectorValues(allVectors_.size());
forAll(allVectors_, vectorI)
{
DynamicList<vectorList>& allTrackVals =
allVectors_[vectorI];
vectorValues[vectorI].setSize(allTrackVals.size());
forAll(allTrackVals, trackI)
{
vectorList& trackVals = allTrackVals[trackI];
vectorValues[vectorI][trackI].transfer(trackVals);
}
}
fileName vtkFile
(
vtkPath
/ vectorFormatterPtr_().getFileName
(
tracks[0],
vectorNames_
)
);
vectorFormatterPtr_().write
(
true, // writeTracks
tracks,
vectorNames_,
vectorValues,
OFstream(vtkFile)()
);
}
}
return true;
}
void Foam::functionObjects::wallBoundedStreamLine::updateMesh
(
const mapPolyMesh& mpm
)
{
if (&mpm.mesh() == &mesh_)
{
read(dict_);
}
}
void Foam::functionObjects::wallBoundedStreamLine::movePoints
(
const polyMesh& mesh
)
{
if (&mesh == &mesh_)
{
// Moving mesh affects the search tree
read(dict_);
}
}
// ************************************************************************* //

View File

@ -1,273 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>.
Class
Foam::functionObjects::wallBoundedStreamLine
Group
grpFieldFunctionObjects
Description
Generates streamline data by sampling a set of user-specified fields along a
particle track, transported by a user-specified velocity field, constrained
to a patch.
Example of function object specification:
\verbatim
wallBoundedStreamLine1
{
type wallBoundedStreamLine;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
setFormat vtk;
U UNear;
trackForward yes;
fields
(
UNear
p
);
lifeTime 10000;
trackLength 1e-3;
nSubCycle 5;
cloudName particleTracks;
seedSampleSet
{
type patchSeed;
patches (wall);
axis x;
maxPoints 20000;
}
}
\endverbatim
Usage
\table
Property | Description | Required | Default value
type | Type name: wallBoundedStreamLine| yes |
setFormat | Output data type | yes |
U | Tracking velocity field name | yes |
fields | Fields to sample | yes |
lifetime | Maximum number of particle tracking steps | yes |
trackLength | Tracking segment length | no |
nSubCycle | Number of tracking steps per cell | no|
cloudName | Cloud name to use | yes |
seedSampleSet| Seeding method (see below)| yes |
\endtable
Where \c seedSampleSet \c type is typically one of
\plaintable
uniform | uniform particle seeding
cloud | cloud of points
patchSeed | seeding via patch faces
triSurfaceMeshPointSet | points according to a tri-surface mesh
\endplaintable
Note
When specifying the track resolution, the \c trackLength OR \c nSubCycle
option should be used
See also
Foam::fvMeshFunctionObject
Foam::functionObjects::timeControl
Foam::sampledSet
Foam::streamLine
SourceFiles
wallBoundedStreamLine.C
\*---------------------------------------------------------------------------*/
#ifndef functionObjects_wallBoundedStreamLine_H
#define functionObjects_wallBoundedStreamLine_H
#include "fvMeshFunctionObject.H"
#include "volFieldsFwd.H"
#include "DynamicList.H"
#include "scalarList.H"
#include "vectorList.H"
#include "writer.H"
#include "indirectPrimitivePatch.H"
#include "tetIndices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class meshSearch;
class sampledSet;
namespace functionObjects
{
/*---------------------------------------------------------------------------*\
Class wallBoundedStreamLine Declaration
\*---------------------------------------------------------------------------*/
class wallBoundedStreamLine
:
public fvMeshFunctionObject
{
// Private data
//- Input dictionary
dictionary dict_;
//- List of fields to sample
wordList fields_;
//- Field to transport particle with
word UName_;
//- Interpolation scheme to use
word interpolationScheme_;
//- Whether to use +u or -u
bool trackForward_;
//- Maximum lifetime (= number of cells) of particle
label lifeTime_;
//- Track length
scalar trackLength_;
//- Optional specified name of particles
word cloudName_;
//- Type of seed
word seedSet_;
//- Names of scalar fields
wordList scalarNames_;
//- Names of vector fields
wordList vectorNames_;
// Demand driven
//- Mesh searching enigne
autoPtr<meshSearch> meshSearchPtr_;
//- Seed set engine
autoPtr<sampledSet> sampledSetPtr_;
//- Axis of the sampled points to output
word sampledSetAxis_;
//- File output writer
autoPtr<writer<scalar>> scalarFormatterPtr_;
autoPtr<writer<vector>> vectorFormatterPtr_;
// Generated data
//- All tracks. Per particle the points it passed through
DynamicList<List<point>> allTracks_;
//- Per scalarField, per particle, the sampled value.
List<DynamicList<scalarList>> allScalars_;
//- Per scalarField, per particle, the sampled value.
List<DynamicList<vectorList>> allVectors_;
//- Construct patch out of all wall patch faces
autoPtr<indirectPrimitivePatch> wallPatch() const;
//- Find wall tet on cell
tetIndices findNearestTet
(
const PackedBoolList& isWallPatch,
const point& seedPt,
const label celli
) const;
//- Do all seeding and tracking
void track();
//- Disallow default bitwise copy construct
wallBoundedStreamLine(const wallBoundedStreamLine&);
//- Disallow default bitwise assignment
void operator=(const wallBoundedStreamLine&);
public:
//- Runtime type information
TypeName("wallBoundedStreamLine");
// Constructors
//- Construct from Time and dictionary
wallBoundedStreamLine
(
const word& name,
const Time& runTime,
const dictionary& dict
);
//- Destructor
virtual ~wallBoundedStreamLine();
// Member Functions
//- Read the field average data
virtual bool read(const dictionary&);
//- Do nothing
virtual bool execute();
//- Calculate and write the wall-bounded streamlines
virtual bool write();
//- Update for changes of mesh
virtual void updateMesh(const mapPolyMesh&);
//- Update for mesh point-motion
virtual void movePoints(const polyMesh&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,435 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "wallBoundedStreamLineParticle.H"
#include "vectorFieldIOField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
(
const trackingData& td,
const point& position,
const label celli,
const label facei
)
{
if (celli == -1)
{
FatalErrorInFunction
<< "Cell:" << celli << abort(FatalError);
}
const tetIndices ti = currentTetIndices();
const vector U = td.vvInterp_[td.UIndex_].interpolate
(
position,
ti, //celli,
facei
);
// Check if at different position
if
(
!sampledPositions_.size()
|| magSqr(sampledPositions_.last()-position) > Foam::sqr(SMALL)
)
{
// Store the location
sampledPositions_.append(position);
// Store the scalar fields
sampledScalars_.setSize(td.vsInterp_.size());
forAll(td.vsInterp_, scalarI)
{
sampledScalars_[scalarI].append
(
td.vsInterp_[scalarI].interpolate
(
position,
ti, //celli,
facei
)
);
}
// Store the vector fields
sampledVectors_.setSize(td.vvInterp_.size());
forAll(td.vvInterp_, vectorI)
{
vector positionU;
if (vectorI == td.UIndex_)
{
positionU = U;
}
else
{
positionU = td.vvInterp_[vectorI].interpolate
(
position,
ti, //celli,
facei
);
}
sampledVectors_[vectorI].append(positionU);
}
}
return U;
}
Foam::vector Foam::wallBoundedStreamLineParticle::sample
(
trackingData& td
)
{
vector U = interpolateFields(td, position(), cell(), tetFace());
if (!td.trackForward_)
{
U = -U;
}
scalar magU = mag(U);
if (magU < SMALL)
{
// Stagnant particle. Might as well stop
lifeTime_ = 0;
}
U /= magU;
return U;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const label tetFacei,
const label tetPtI,
const label meshEdgeStart,
const label diagEdge,
const label lifeTime
)
:
wallBoundedParticle
(
mesh,
position,
celli,
tetFacei,
tetPtI,
meshEdgeStart,
diagEdge
),
lifeTime_(lifeTime)
{}
Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
)
:
wallBoundedParticle(mesh, is, readFields)
{
if (readFields)
{
List<scalarList> sampledScalars;
List<vectorList> sampledVectors;
is >> lifeTime_
>> sampledPositions_ >> sampledScalars >> sampledVectors;
sampledScalars_.setSize(sampledScalars.size());
forAll(sampledScalars, i)
{
sampledScalars_[i].transfer(sampledScalars[i]);
}
sampledVectors_.setSize(sampledVectors.size());
forAll(sampledVectors, i)
{
sampledVectors_[i].transfer(sampledVectors[i]);
}
}
// Check state of Istream
is.check
(
"wallBoundedStreamLineParticle::wallBoundedStreamLineParticle"
"(const Cloud<wallBoundedStreamLineParticle>&, Istream&, bool)"
);
}
Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
(
const wallBoundedStreamLineParticle& p
)
:
wallBoundedParticle(p),
lifeTime_(p.lifeTime_),
sampledPositions_(p.sampledPositions_),
sampledScalars_(p.sampledScalars_),
sampledVectors_(p.sampledVectors_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::wallBoundedStreamLineParticle::move
(
trackingData& td,
const scalar trackTime
)
{
wallBoundedStreamLineParticle& p = static_cast
<
wallBoundedStreamLineParticle&
>(*this);
// Check position is inside tet
//checkInside();
td.switchProcessor = false;
td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar maxDt = mesh_.bounds().mag();
while
(
td.keepParticle
&& !td.switchProcessor
&& lifeTime_ > 0
)
{
// set the lagrangian time-step
scalar dt = maxDt;
--lifeTime_;
// Get sampled velocity and fields. Store if position changed.
vector U = sample(td);
// !user parameter!
if (dt < SMALL)
{
// Force removal
lifeTime_ = 0;
break;
}
if (td.trackLength_ < GREAT)
{
dt = td.trackLength_;
}
scalar fraction = trackToEdge(td, position() + dt*U);
dt *= fraction;
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
if (tEnd <= ROOTVSMALL)
{
// Force removal
lifeTime_ = 0;
}
if
(
!td.keepParticle
|| td.switchProcessor
|| lifeTime_ == 0
)
{
break;
}
}
if (!td.keepParticle || lifeTime_ == 0)
{
if (lifeTime_ == 0)
{
if (debug)
{
Pout<< "wallBoundedStreamLineParticle :"
<< " Removing stagnant particle:"
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
td.keepParticle = false;
}
else
{
// Normal exit. Store last position and fields
sample(td);
if (debug)
{
Pout<< "wallBoundedStreamLineParticle : Removing particle:"
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
}
// Transfer particle data into trackingData.
{
//td.allPositions_.append(sampledPositions_);
td.allPositions_.append(vectorList());
vectorList& top = td.allPositions_.last();
top.transfer(sampledPositions_);
}
forAll(sampledScalars_, i)
{
//td.allScalars_[i].append(sampledScalars_[i]);
td.allScalars_[i].append(scalarList());
scalarList& top = td.allScalars_[i].last();
top.transfer(sampledScalars_[i]);
}
forAll(sampledVectors_, i)
{
//td.allVectors_[i].append(sampledVectors_[i]);
td.allVectors_[i].append(vectorList());
vectorList& top = td.allVectors_[i].last();
top.transfer(sampledVectors_[i]);
}
}
return td.keepParticle;
}
void Foam::wallBoundedStreamLineParticle::readFields
(
Cloud<wallBoundedStreamLineParticle>& c
)
{
if (!c.size())
{
return;
}
wallBoundedParticle::readFields(c);
IOField<label> lifeTime
(
c.fieldIOobject("lifeTime", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, lifeTime);
vectorFieldIOField sampledPositions
(
c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, sampledPositions);
label i = 0;
forAllIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
{
iter().lifeTime_ = lifeTime[i];
iter().sampledPositions_.transfer(sampledPositions[i]);
i++;
}
}
void Foam::wallBoundedStreamLineParticle::writeFields
(
const Cloud<wallBoundedStreamLineParticle>& c
)
{
wallBoundedParticle::writeFields(c);
label np = c.size();
IOField<label> lifeTime
(
c.fieldIOobject("lifeTime", IOobject::NO_READ),
np
);
vectorFieldIOField sampledPositions
(
c.fieldIOobject("sampledPositions", IOobject::NO_READ),
np
);
label i = 0;
forAllConstIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
{
lifeTime[i] = iter().lifeTime_;
sampledPositions[i] = iter().sampledPositions_;
i++;
}
lifeTime.write();
sampledPositions.write();
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const wallBoundedStreamLineParticle& p
)
{
os << static_cast<const wallBoundedParticle&>(p)
<< token::SPACE << p.lifeTime_
<< token::SPACE << p.sampledPositions_
<< token::SPACE << p.sampledScalars_
<< token::SPACE << p.sampledVectors_;
// Check state of Ostream
os.check
(
"Ostream& operator<<(Ostream&, const wallBoundedStreamLineParticle&)"
);
return os;
}
// ************************************************************************* //

View File

@ -1,260 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
Class
Foam::wallBoundedStreamLineParticle
Description
Particle class that samples fields as it passes through. Used in streamline
calculation.
SourceFiles
wallBoundedStreamLineParticle.C
\*---------------------------------------------------------------------------*/
#ifndef wallBoundedStreamLineParticle_H
#define wallBoundedStreamLineParticle_H
#include "wallBoundedParticle.H"
#include "autoPtr.H"
#include "interpolation.H"
#include "vectorList.H"
#include "InfoProxy.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class wallBoundedStreamLineParticleCloud;
// Forward declaration of friend functions and operators
class wallBoundedStreamLineParticle;
Ostream& operator<<(Ostream&, const wallBoundedStreamLineParticle&);
/*---------------------------------------------------------------------------*\
Class wallBoundedStreamLineParticle Declaration
\*---------------------------------------------------------------------------*/
class wallBoundedStreamLineParticle
:
public wallBoundedParticle
{
public:
//- Class used to pass tracking data to the trackToEdge function
class trackingData
:
public wallBoundedParticle::TrackingData
<
Cloud<wallBoundedStreamLineParticle>
>
{
public:
const PtrList<interpolation<scalar>>& vsInterp_;
const PtrList<interpolation<vector>>& vvInterp_;
const label UIndex_;
const bool trackForward_;
const scalar trackLength_;
DynamicList<vectorList>& allPositions_;
List<DynamicList<scalarList>>& allScalars_;
List<DynamicList<vectorList>>& allVectors_;
// Constructors
trackingData
(
Cloud<wallBoundedStreamLineParticle>& cloud,
const PtrList<interpolation<scalar>>& vsInterp,
const PtrList<interpolation<vector>>& vvInterp,
const label UIndex,
const bool trackForward,
const scalar trackLength,
const PackedBoolList& isWallPatch,
DynamicList<List<point>>& allPositions,
List<DynamicList<scalarList>>& allScalars,
List<DynamicList<vectorList>>& allVectors
)
:
wallBoundedParticle::TrackingData
<
Cloud<wallBoundedStreamLineParticle>
>
(
cloud,
isWallPatch
),
vsInterp_(vsInterp),
vvInterp_(vvInterp),
UIndex_(UIndex),
trackForward_(trackForward),
trackLength_(trackLength),
allPositions_(allPositions),
allScalars_(allScalars),
allVectors_(allVectors)
{}
};
private:
// Private data
//- Lifetime of particle. Particle dies when reaches 0.
label lifeTime_;
//- Sampled positions
DynamicList<point> sampledPositions_;
//- Sampled scalars
List<DynamicList<scalar>> sampledScalars_;
//- Sampled vectors
List<DynamicList<vector>> sampledVectors_;
// Private Member Functions
vector interpolateFields
(
const trackingData& td,
const point& position,
const label celli,
const label facei
);
vector sample(trackingData& td);
public:
// Constructors
//- Construct from components
wallBoundedStreamLineParticle
(
const polyMesh& c,
const vector& position,
const label celli,
const label tetFacei,
const label tetPtI,
const label meshEdgeStart,
const label diagEdge,
const label lifeTime
);
//- Construct from Istream
wallBoundedStreamLineParticle
(
const polyMesh& c,
Istream& is,
bool readFields = true
);
//- Construct copy
wallBoundedStreamLineParticle(const wallBoundedStreamLineParticle& p);
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new wallBoundedStreamLineParticle(*this));
}
//- Factory class to read-construct particles used for
// parallel transfer
class iNew
{
const polyMesh& mesh_;
public:
iNew(const polyMesh& mesh)
:
mesh_(mesh)
{}
autoPtr<wallBoundedStreamLineParticle> operator()
(
Istream& is
) const
{
return autoPtr<wallBoundedStreamLineParticle>
(
new wallBoundedStreamLineParticle(mesh_, is, true)
);
}
};
// Member Functions
// Tracking
//- Track all particles to their end point
bool move(trackingData&, const scalar trackTime);
// I-O
//- Read
static void readFields(Cloud<wallBoundedStreamLineParticle>&);
//- Write
static void writeFields
(
const Cloud<wallBoundedStreamLineParticle>&
);
// Ostream Operator
friend Ostream& operator<<
(
Ostream&,
const wallBoundedStreamLineParticle&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -176,16 +176,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
U += velocity; U += velocity;
addNewParcel addNewParcel(p, celli, U, Ei, typeId);
(
p,
U,
Ei,
celli,
cellTetIs.face(),
cellTetIs.tetPt(),
typeId
);
} }
} }
} }
@ -464,27 +455,13 @@ template<class ParcelType>
void Foam::DSMCCloud<ParcelType>::addNewParcel void Foam::DSMCCloud<ParcelType>::addNewParcel
( (
const vector& position, const vector& position,
const label celli,
const vector& U, const vector& U,
const scalar Ei, const scalar Ei,
const label celli,
const label tetFacei,
const label tetPti,
const label typeId const label typeId
) )
{ {
ParcelType* pPtr = new ParcelType this->addParticle(new ParcelType(mesh_, position, celli, U, Ei, typeId));
(
mesh_,
position,
U,
Ei,
celli,
tetFacei,
tetPti,
typeId
);
this->addParticle(pPtr);
} }
@ -1112,9 +1089,7 @@ void Foam::DSMCCloud<ParcelType>::dumpParticlePositions() const
template<class ParcelType> template<class ParcelType>
void Foam::DSMCCloud<ParcelType>::autoMap(const mapPolyMesh& mapper) void Foam::DSMCCloud<ParcelType>::autoMap(const mapPolyMesh& mapper)
{ {
typedef typename ParcelType::trackingData tdType; Cloud<ParcelType>::autoMap(mapper);
tdType td(*this);
Cloud<ParcelType>::template autoMap<tdType>(td, mapper);
// Update the cell occupancy field // Update the cell occupancy field
cellOccupancy_.setSize(mesh_.nCells()); cellOccupancy_.setSize(mesh_.nCells());

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -44,6 +44,7 @@ SourceFiles
#include "fvMesh.H" #include "fvMesh.H"
#include "volFields.H" #include "volFields.H"
#include "scalarIOField.H" #include "scalarIOField.H"
#include "barycentric.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -451,11 +452,9 @@ public:
void addNewParcel void addNewParcel
( (
const vector& position, const vector& position,
const label celli,
const vector& U, const vector& U,
const scalar Ei, const scalar Ei,
const label celli,
const label tetFacei,
const label tetPti,
const label typeId const label typeId
); );

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -41,19 +41,16 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
const polyMesh& mesh = td.cloud().pMesh(); const polyMesh& mesh = td.cloud().pMesh();
const polyBoundaryMesh& pbMesh = mesh.boundaryMesh(); const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
scalar tEnd = (1.0 - p.stepFraction())*trackTime;
const scalar dtMax = tEnd;
// For reduced-D cases, the velocity used to track needs to be // For reduced-D cases, the velocity used to track needs to be
// constrained, but the actual U_ of the parcel must not be // constrained, but the actual U_ of the parcel must not be
// altered or used, as it is altered by patch interactions an // altered or used, as it is altered by patch interactions an
// needs to retain its 3D value for collision purposes. // needs to retain its 3D value for collision purposes.
vector Utracking = U_; vector Utracking = U_;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL) while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
{ {
// Apply correction to position for reduced-D cases // Apply correction to position for reduced-D cases
meshTools::constrainToMeshCentre(mesh, p.position()); p.constrainToMeshCentre();
Utracking = U_; Utracking = U_;
@ -61,16 +58,10 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
// reduced-D cases // reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking); meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
// Set the Lagrangian time-step const scalar f = 1 - p.stepFraction();
scalar dt = min(dtMax, tEnd); p.trackToFace(f*trackTime*Utracking, f, td);
dt *= p.trackToFace(p.position() + dt*Utracking, td); if (p.onBoundaryFace() && td.keepParticle)
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/trackTime;
if (p.onBoundary() && td.keepParticle)
{ {
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())])) if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{ {

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -178,12 +178,24 @@ public:
inline DSMCParcel inline DSMCParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const vector& U,
const scalar Ei,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
const vector& U,
const scalar Ei,
const label typeId
);
//- Construct from a position and a cell, searching for the rest of the
// required topology
inline DSMCParcel
(
const polyMesh& mesh,
const vector& position,
const label celli,
const vector& U,
const scalar Ei,
const label typeId const label typeId
); );

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -55,16 +55,34 @@ template<class ParcelType>
inline Foam::DSMCParcel<ParcelType>::DSMCParcel inline Foam::DSMCParcel<ParcelType>::DSMCParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const vector& U,
const scalar Ei,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
const vector& U,
const scalar Ei,
const label typeId const label typeId
) )
: :
ParcelType(mesh, position, celli, tetFacei, tetPti), ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
U_(U),
Ei_(Ei),
typeId_(typeId)
{}
template<class ParcelType>
inline Foam::DSMCParcel<ParcelType>::DSMCParcel
(
const polyMesh& mesh,
const vector& position,
const label celli,
const vector& U,
const scalar Ei,
const label typeId
)
:
ParcelType(mesh, position, celli),
U_(U), U_(U),
Ei_(Ei), Ei_(Ei),
typeId_(typeId) typeId_(typeId)

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -403,16 +403,7 @@ void Foam::FreeStream<CloudType>::inflow()
cloud.constProps(typeId).internalDegreesOfFreedom() cloud.constProps(typeId).internalDegreesOfFreedom()
); );
cloud.addNewParcel cloud.addNewParcel(p, celli, U, Ei, typeId);
(
p,
U,
Ei,
celli,
globalFaceIndex,
faceTetIs.tetPt(),
typeId
);
particlesInserted++; particlesInserted++;
} }

View File

@ -105,7 +105,8 @@ Foam::Cloud<ParticleType>::Cloud
polyMesh_(pMesh), polyMesh_(pMesh),
labels_(), labels_(),
nTrackingRescues_(), nTrackingRescues_(),
cellWallFacesPtr_() cellWallFacesPtr_(),
globalPositionsPtr_()
{ {
checkPatches(); checkPatches();
@ -234,6 +235,8 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
// Allocate transfer buffers // Allocate transfer buffers
PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking); PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
// Clear the global positions as there are about to change
globalPositionsPtr_.clear();
// While there are particles to transfer // While there are particles to transfer
while (true) while (true)
@ -393,21 +396,16 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
template<class ParticleType> template<class ParticleType>
template<class TrackData> void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
void Foam::Cloud<ParticleType>::autoMap
(
TrackData& td,
const mapPolyMesh& mapper
)
{ {
if (cloud::debug) if (!globalPositionsPtr_.valid())
{ {
InfoInFunction << "for lagrangian cloud " << cloud::name() << endl; FatalErrorInFunction
<< "Global positions are not available. "
<< "Cloud::storeGlobalPositions has not been called."
<< exit(FatalError);
} }
const labelList& reverseCellMap = mapper.reverseCellMap();
const labelList& reverseFaceMap = mapper.reverseFaceMap();
// Reset stored data that relies on the mesh // Reset stored data that relies on the mesh
// polyMesh_.clearCellTree(); // polyMesh_.clearCellTree();
cellWallFacesPtr_.clear(); cellWallFacesPtr_.clear();
@ -417,51 +415,13 @@ void Foam::Cloud<ParticleType>::autoMap
// there is a comms mismatch. // there is a comms mismatch.
polyMesh_.tetBasePtIs(); polyMesh_.tetBasePtIs();
const vectorField& positions = globalPositionsPtr_();
forAllIter(typename Cloud<ParticleType>, *this, pIter) label i = 0;
forAllIter(typename Cloud<ParticleType>, *this, iter)
{ {
ParticleType& p = pIter(); iter().autoMap(positions[i], mapper);
++ i;
if (reverseCellMap[p.cell()] >= 0)
{
p.cell() = reverseCellMap[p.cell()];
if (p.face() >= 0 && reverseFaceMap[p.face()] >= 0)
{
p.face() = reverseFaceMap[p.face()];
}
else
{
p.face() = -1;
}
p.initCellFacePt();
}
else
{
label trackStartCell = mapper.mergedCell(p.cell());
if (trackStartCell < 0)
{
trackStartCell = 0;
p.cell() = 0;
}
else
{
p.cell() = trackStartCell;
}
vector pos = p.position();
const_cast<vector&>(p.position()) =
polyMesh_.cellCentres()[trackStartCell];
p.stepFraction() = 0;
p.initCellFacePt();
p.track(pos, td);
}
} }
} }
@ -485,6 +445,27 @@ void Foam::Cloud<ParticleType>::writePositions() const
} }
template<class ParticleType>
void Foam::Cloud<ParticleType>::storeGlobalPositions() const
{
// Store the global positions for later use by autoMap. It would be
// preferable not to need this. If the mapPolyMesh object passed to autoMap
// had a copy of the old mesh then the global positions could be recovered
// within autoMap, and this pre-processing would not be necessary.
globalPositionsPtr_.reset(new vectorField(this->size()));
vectorField& positions = globalPositionsPtr_();
label i = 0;
forAllConstIter(typename Cloud<ParticleType>, *this, iter)
{
positions[i] = iter().position();
++ i;
}
}
// * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * // // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
#include "CloudIO.C" #include "CloudIO.C"

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -87,6 +87,9 @@ class Cloud
//- Does the cell have wall faces //- Does the cell have wall faces
mutable autoPtr<PackedBoolList> cellWallFacesPtr_; mutable autoPtr<PackedBoolList> cellWallFacesPtr_;
//- Temporary storage for the global particle positions
mutable autoPtr<vectorField> globalPositionsPtr_;
// Private Member Functions // Private Member Functions
@ -162,7 +165,7 @@ public:
return IDLList<ParticleType>::size(); return IDLList<ParticleType>::size();
}; };
DynamicList<label>& labels() DynamicList<label>& labels() const
{ {
return labels_; return labels_;
} }
@ -256,8 +259,7 @@ public:
//- Remap the cells of particles corresponding to the //- Remap the cells of particles corresponding to the
// mesh topology change // mesh topology change
template<class TrackData> void autoMap(const mapPolyMesh&);
void autoMap(TrackData& td, const mapPolyMesh&);
// Read // Read
@ -304,6 +306,10 @@ public:
//- Write positions to \<cloudName\>_positions.obj file //- Write positions to \<cloudName\>_positions.obj file
void writePositions() const; void writePositions() const;
//- Call this before a topology change. Stores the particles global
// positions in the database for use during mapping.
void storeGlobalPositions() const;
// Ostream Operator // Ostream Operator

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -132,13 +132,6 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
// them, otherwise, if some processors have no particles then // them, otherwise, if some processors have no particles then
// there is a comms mismatch. // there is a comms mismatch.
polyMesh_.tetBasePtIs(); polyMesh_.tetBasePtIs();
forAllIter(typename Cloud<ParticleType>, *this, pIter)
{
ParticleType& p = pIter();
p.initCellFacePt();
}
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -965,14 +965,7 @@ void Foam::InteractionLists<ParticleType>::prepareParticleToBeReferred
globalTransforms.transformIndex(ciat) globalTransforms.transformIndex(ciat)
); );
particle->position() = transform.invTransformPosition(particle->position()); particle->prepareForInteractionListReferral(transform);
particle->transformProperties(-transform.t());
if (transform.hasR())
{
particle->transformProperties(transform.R().T());
}
} }
@ -1230,6 +1223,15 @@ void Foam::InteractionLists<ParticleType>::receiveReferredData
} }
} }
forAll(referredParticles_, refCelli)
{
IDLList<ParticleType>& refCell = referredParticles_[refCelli];
forAllIter(typename IDLList<ParticleType>, refCell, iter)
{
iter().correctAfterInteractionListReferral(ril_[refCelli][0]);
}
}
fillReferredParticleCloud(); fillReferredParticleCloud();
wallFaceMap().receive(pBufs, referredWallData_); wallFaceMap().receive(pBufs, referredWallData_);

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -65,27 +65,14 @@ public:
indexedParticle indexedParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
const label index = 0 const label index = 0
) )
: :
particle(mesh, position, celli, tetFacei, tetPti), particle(mesh, coordinates, celli, tetFacei, tetPti),
index_(index)
{}
//- Construct from components, with searching for tetFace and tetPt
indexedParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const label index = 0
)
:
particle(mesh, position, celli),
index_(index) index_(index)
{} {}

File diff suppressed because it is too large Load Diff

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -33,6 +33,8 @@ Description
#define particle_H #define particle_H
#include "vector.H" #include "vector.H"
#include "barycentric.H"
#include "barycentricTensor.H"
#include "Cloud.H" #include "Cloud.H"
#include "IDLList.H" #include "IDLList.H"
#include "pointField.H" #include "pointField.H"
@ -42,6 +44,7 @@ Description
#include "FixedList.H" #include "FixedList.H"
#include "polyMeshTetDecomposition.H" #include "polyMeshTetDecomposition.H"
#include "particleMacros.H" #include "particleMacros.H"
#include "vectorTensorTransform.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -130,25 +133,19 @@ public:
}; };
protected: private:
// Protected data // Private data
//- Reference to the polyMesh database //- Reference to the polyMesh database
const polyMesh& mesh_; const polyMesh& mesh_;
//- Position of particle //- Coordinates of particle
vector position_; barycentric coordinates_;
//- Index of the cell it is in //- Index of the cell it is in
label celli_; label celli_;
//- Face index if the particle is on a face otherwise -1
label facei_;
//- Fraction of time-step completed
scalar stepFraction_;
//- Index of the face that owns the decomposed tet that the //- Index of the face that owns the decomposed tet that the
// particle is in // particle is in
label tetFacei_; label tetFacei_;
@ -158,6 +155,12 @@ protected:
// point. // point.
label tetPti_; label tetPti_;
//- Face index if the particle is on a face otherwise -1
label facei_;
//- Fraction of time-step completed
scalar stepFraction_;
//- Originating processor id //- Originating processor id
label origProc_; label origProc_;
@ -165,74 +168,137 @@ protected:
label origId_; label origId_;
private:
// Private Member Functions // Private Member Functions
//- Find the tet tri faces between position and tet centre // Tetrahedra functions
void findTris
(
const vector& position,
DynamicList<label>& faceList,
const tetPointRef& tet,
const FixedList<vector, 4>& tetAreas,
const FixedList<label, 4>& tetPlaneBasePtIs,
const scalar tol
) const;
//- Find the lambda value for the line to-from across the //- Get indices into the current face for the face-bound vertices of
// given tri face, where p = from + lambda*(to - from) // the current tet.
inline scalar tetLambda void tetFaceIndices
( (
const vector& from, label& baseI,
const vector& to, label& vertex1I,
const label triI, label& vertex2I
const vector& tetArea, ) const;
const label tetPlaneBasePtI,
const label celli,
const label tetFacei,
const label tetPti,
const scalar tol
) const;
//- Find the lambda value for a moving tri face //- Get indices into the mesh points for the face-bound vertices of
inline scalar movingTetLambda // the current tet.
( void tetMeshIndices
const vector& from, (
const vector& to, label& baseI,
const label triI, label& vertex1I,
const vector& tetArea, label& vertex2I
const label tetPlaneBasePtI, ) const;
const label celli,
const label tetFacei,
const label tetPti,
const scalar tol
) const;
//- Modify the tet owner data by crossing triI //- Get the vertices of the current tet
inline void tetNeighbour(label triI); void tetGeometry
(
vector& centre,
vector& base,
vector& vertex1,
vector& vertex2
) const;
//- Cross the from the given face across the given edge of the //- Get the transformation associated with the current tet. This
// given cell to find the resulting face and tetPti // will convert a barycentric position within the tet to a
inline void crossEdgeConnectedFace // cartesian position in the global coordinate system. The
( // conversion is x = A & y, where x is the cartesian position, y is
const label& celli, // the barycentric position and A is the transformation tensor.
label& tetFacei, barycentricTensor tetTransform() const;
label& tetPti,
const edge& e
);
//- Hit wall faces in the current cell if the //- Get the reverse transform associated with the current tet. The
//- wallImpactDistance is non-zero. They may not be in // conversion is detA*y = (x - centre) & T. The variables x, y and
//- Different tets to the current. // centre have the same meaning as for the forward transform. T is
template<class CloudType> // the transposed inverse of the forward transform tensor, A,
inline void hitWallFaces // multiplied by its determinant, detA. This separation allows
( // the barycentric tracking algorithm to function on inverted or
const CloudType& td, // degenerate tetrahedra.
const vector& from, void tetReverseTransform
const vector& to, (
scalar& lambdaMin, vector& centre,
tetIndices& closestTetIs scalar& detA,
); barycentricTensor& T
) const;
//- 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
(
const scalar endStepFraction,
Pair<vector>& centre,
Pair<vector>& base,
Pair<vector>& vertex1,
Pair<vector>& vertex2
) const;
//- Get the transformation associated with the current, moving, tet.
// 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<barycentricTensor> movingTetTransform
(
const scalar endStepFraction
) const;
//- Get the reverse transformation associated with the current,
// moving, tet. This is of the same form as for the static case. As
// with the moving geometry, a function of the tracking fraction is
// returned for each component. The functions are higher order than
// for the forward transform; the determinant is cubic, and the
// tensor is quadratic.
void movingTetReverseTransform
(
const scalar endStepFraction,
Pair<vector>& centre,
FixedList<scalar, 4>& detA,
FixedList<barycentricTensor, 3>& T
) const;
// Transformations
//- Reflection transform. Corrects the coordinates when the particle
// moves between two tets which share a base vertex, but for which
// the other two non cell-centre vertices are reversed. All hits
// which retain the same face behave this way, as do face hits.
void reflect();
//- Rotation transform. Corrects the coordinates when the particle
// moves between two tets with different base vertices, but are
// otherwise similarly oriented. Hits which change the face within
// the cell make use of both this and the reflect transform.
void rotate(const bool direction);
// Topology changes
//- Change tet within a cell. Called after a triangle is hit.
void changeTet(const label tetTriI);
//- Change tet face within a cell. Called by changeTet.
void changeFace(const label tetTriI);
//- Change cell. Called when the particle hits an internal face.
void changeCell();
// Geometry changes
//- Locate the particle at the given position
void locate
(
const vector& position,
const vector* direction,
const label celli,
const bool boundaryFail,
const string boundaryMsg
);
protected:
// Patch interactions // Patch interactions
@ -315,24 +381,13 @@ public:
//- String representation of properties //- String representation of properties
DefinePropertyList DefinePropertyList
( (
"(Px Py Pz) celli facei stepFraction " "(Px Py Pz) celli tetFacei tetPti "
"tetFacei tetPti origProc origId" "facei stepFraction origProc origId"
); );
//- Cumulative particle counter - used to provode unique ID //- Cumulative particle counter - used to provode unique ID
static label particleCount_; static label particleCount_;
//- Fraction of distance to tet centre to move a particle to
// 'rescue' it from a tracking problem
static const scalar trackingCorrectionTol;
//- Fraction of the cell volume to use in determining tolerance values
// for the denominator and numerator of lambda
static const scalar lambdaDistanceToleranceCoeff;
//- Minimum stepFraction tolerance
static const scalar minStepFractionTol;
// Constructors // Constructors
@ -340,20 +395,19 @@ public:
particle particle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from components, tetFacei_ and tetPti_ are not //- Construct from a position and a cell, searching for the rest of the
// supplied so they will be deduced by a search // required topology
particle particle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli, const label celli
bool doCellFacePt = true
); );
//- Construct from Istream //- Construct from Istream
@ -406,27 +460,15 @@ public:
//- Return the mesh database //- Return the mesh database
inline const polyMesh& mesh() const; inline const polyMesh& mesh() const;
//- Return current particle position //- Return current particle coordinates
inline const vector& position() const; inline const barycentric& coordinates() const;
//- Return current particle position
inline vector& position();
//- Return current cell particle is in
inline label& cell();
//- Return current cell particle is in //- Return current cell particle is in
inline label cell() const; inline label cell() const;
//- Return current tet face particle is in
inline label& tetFace();
//- Return current tet face particle is in //- Return current tet face particle is in
inline label tetFace() const; inline label tetFace() const;
//- Return current tet face particle is in
inline label& tetPt();
//- Return current tet face particle is in //- Return current tet face particle is in
inline label tetPt() const; inline label tetPt() const;
@ -447,43 +489,32 @@ public:
// on oldPoints // on oldPoints
inline vector oldNormal() const; inline vector oldNormal() const;
//- Return current face particle is on otherwise -1
inline label& face();
//- Return current face particle is on otherwise -1 //- Return current face particle is on otherwise -1
inline label face() const; inline label face() const;
//- Return the impact model to be used, soft or hard (default).
inline bool softImpact() const;
//- Return the particle current time //- Return the particle current time
inline scalar currentTime() const; inline scalar currentTime() const;
// Check // Check
//- Check the stored cell value (setting if necessary) and //- Is the particle on a face?
// initialise the tetFace and tetPt values inline bool onFace() const;
inline void initCellFacePt();
//- Is the particle on the boundary/(or outside the domain)? //- Is the particle on an internal face?
inline bool onBoundary() const; inline bool onInternalFace() const;
//- Is this global face an internal face? //- Is the particle on a boundary face?
inline bool internalFace(const label facei) const; inline bool onBoundaryFace() const;
//- Is this global face a boundary face? //- Which patch is particle on // <-- !!!
inline bool boundaryFace(const label facei) const;
//- Which patch is particle on
inline label patch(const label facei) const; inline label patch(const label facei) const;
//- Which face of this patch is this particle on //- Which face of this patch is this particle on // <-- !!!
inline label patchFace inline label patchFace(const label patchi, const label facei) const;
(
const label patchi, //- Return current particle position
const label facei inline vector position() const;
) const;
//- Return the fraction of time-step completed //- Return the fraction of time-step completed
inline scalar& stepFraction(); inline scalar& stepFraction();
@ -506,28 +537,57 @@ public:
// Track // Track
//- Track particle to end of trajectory //- Track along the displacement for a given fraction of the overall
// or until it hits the boundary. // step. End when the track is complete, or when a boundary is hit.
// On entry 'stepFraction()' should be set to the fraction of the // On exit, stepFraction_ will have been incremented to the current
// time-step at which the tracking starts and on exit it contains // position, and facei_ will be set to the index of the boundary
// the fraction of the time-step completed. // face that was hit, or -1 if the track completed within a cell.
// Returns the boundary face index if the track stops at the // The proportion of the displacement still to be completed is
// boundary, -1 otherwise. // returned.
template<class TrackData> scalar track
label track(const vector& endPosition, TrackData& td); (
const vector& displacement,
const scalar fraction
);
//- Track particle to a given position and returns 1.0 if the //- As particle::track, but also stops on internal faces.
// trajectory is completed without hitting a face otherwise scalar trackToFace
// stops at the face and returns the fraction of the trajectory (
// completed. const vector& displacement,
// on entry 'stepFraction()' should be set to the fraction of the const scalar fraction
// time-step at which the tracking starts. );
template<class TrackData>
scalar trackToFace(const vector& endPosition, TrackData& td);
//- Return the index of the face to be used in the interpolation //- As particle::trackToFace, but also stops on tet triangles. On
// routine // exit, tetTriI is set to the index of the tet triangle that was
inline label faceInterpolation() const; // hit, or -1 if the end position was reached.
scalar trackToStationaryTri
(
const vector& displacement,
const scalar fraction,
label& tetTriI
);
//- As particle::trackToTri, but for moving meshes
scalar trackToMovingTri
(
const vector& displacement,
const scalar fraction,
label& tetTriI
);
//- As non-templated particle::trackToFace, but with additional
// boundary handling.
template<class TrackData>
void trackToFace
(
const vector& displacement,
const scalar fraction,
TrackData& td
);
//- Set the constrained components of the particle position to the
// mesh centre.
void constrainToMeshCentre();
// Transformations // Transformations
@ -558,6 +618,38 @@ public:
void correctAfterParallelTransfer(const label patchi, TrackData& td); void correctAfterParallelTransfer(const label patchi, TrackData& td);
// Interaction list referral
//- Break the topology and store the particle position so that the
// particle can be referred.
void prepareForInteractionListReferral
(
const vectorTensorTransform& transform
);
//- Correct the topology after referral. The particle may still be
// outside the stored tet and therefore not track-able.
void correctAfterInteractionListReferral(const label celli);
// Decompose and reconstruct
//- Return the tet point approproate for decomposition or reconstruction
// to or from the given mesh.
label procTetPt
(
const polyMesh& procMesh,
const label procCell,
const label procTetFace
) const;
// Mapping
//- Map after a topology change
void autoMap(const vector& position, const mapPolyMesh& mapper);
// I-O // I-O
//- Read the fields associated with the owner cloud //- Read the fields associated with the owner cloud

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -26,541 +26,6 @@ License
#include "polyMesh.H" #include "polyMesh.H"
#include "Time.H" #include "Time.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::particle::findTris
(
const vector& position,
DynamicList<label>& faceList,
const tetPointRef& tet,
const FixedList<vector, 4>& tetAreas,
const FixedList<label, 4>& tetPlaneBasePtIs,
const scalar tol
) const
{
faceList.clear();
const point Ct = tet.centre();
for (label i = 0; i < 4; i++)
{
scalar lambda = tetLambda
(
Ct,
position,
i,
tetAreas[i],
tetPlaneBasePtIs[i],
celli_,
tetFacei_,
tetPti_,
tol
);
if ((lambda > 0.0) && (lambda < 1.0))
{
faceList.append(i);
}
}
}
inline Foam::scalar Foam::particle::tetLambda
(
const vector& from,
const vector& to,
const label triI,
const vector& n,
const label tetPlaneBasePtI,
const label celli,
const label tetFacei,
const label tetPti,
const scalar tol
) const
{
const pointField& pPts = mesh_.points();
if (mesh_.moving())
{
return movingTetLambda
(
from,
to,
triI,
n,
tetPlaneBasePtI,
celli,
tetFacei,
tetPti,
tol
);
}
const point& base = pPts[tetPlaneBasePtI];
scalar lambdaNumerator = (base - from) & n;
scalar lambdaDenominator = (to - from) & n;
// n carries the area of the tet faces, so the dot product with a
// delta-length has the units of volume. Comparing the component of each
// delta-length in the direction of n times the face area to a fraction of
// the cell volume.
if (mag(lambdaDenominator) < tol)
{
if (mag(lambdaNumerator) < tol)
{
// Track starts on the face, and is potentially
// parallel to it. +-tol/+-tol is not a good
// comparison, return 0.0, in anticipation of tet
// centre correction.
return 0.0;
}
else
{
if (mag((to - from)) < tol/mag(n))
{
// 'Zero' length track (compared to the tolerance, which is
// based on the cell volume, divided by the tet face area), not
// along the face, face cannot be crossed.
return GREAT;
}
else
{
// Trajectory is non-zero and parallel to face
lambdaDenominator = sign(lambdaDenominator)*SMALL;
}
}
}
return lambdaNumerator/lambdaDenominator;
}
inline Foam::scalar Foam::particle::movingTetLambda
(
const vector& from,
const vector& to,
const label triI,
const vector& n,
const label tetPlaneBasePtI,
const label celli,
const label tetFacei,
const label tetPti,
const scalar tol
) const
{
const pointField& pPts = mesh_.points();
const pointField& oldPPts = mesh_.oldPoints();
// Base point of plane at end of motion
const point& b = pPts[tetPlaneBasePtI];
// n: Normal of plane at end of motion
// Base point of plane at start of timestep
const point& b00 = oldPPts[tetPlaneBasePtI];
// Base point of plane at start of tracking portion (cast forward by
// stepFraction)
point b0 = b00 + stepFraction_*(b - b00);
// Normal of plane at start of tracking portion
vector n0 = Zero;
{
tetIndices tetIs(celli, tetFacei, tetPti, mesh_);
// Tet at timestep start
tetPointRef tet00 = tetIs.oldTet(mesh_);
// Tet at timestep end
tetPointRef tet = tetIs.tet(mesh_);
point tet0PtA = tet00.a() + stepFraction_*(tet.a() - tet00.a());
point tet0PtB = tet00.b() + stepFraction_*(tet.b() - tet00.b());
point tet0PtC = tet00.c() + stepFraction_*(tet.c() - tet00.c());
point tet0PtD = tet00.d() + stepFraction_*(tet.d() - tet00.d());
// Tracking portion start tet (cast forward by stepFraction)
tetPointRef tet0(tet0PtA, tet0PtB, tet0PtC, tet0PtD);
switch (triI)
{
case 0:
{
n0 = tet0.Sa();
break;
}
case 1:
{
n0 = tet0.Sb();
break;
}
case 2:
{
n0 = tet0.Sc();
break;
}
case 3:
{
n0 = tet0.Sd();
break;
}
default:
{
break;
}
}
}
if (mag(n0) < SMALL)
{
// If the old normal is zero (for example in layer addition)
// then use the current normal;
n0 = n;
}
scalar lambdaNumerator = 0;
scalar lambdaDenominator = 0;
vector dP = to - from;
vector dN = n - n0;
vector dB = b - b0;
vector dS = from - b0;
if (mag(dN) > SMALL)
{
scalar a = (dP - dB) & dN;
scalar b = ((dP - dB) & n0) + (dS & dN);
scalar c = dS & n0;
if (mag(a) > SMALL)
{
// Solve quadratic for lambda
scalar discriminant = sqr(b) - 4.0*a*c;
if (discriminant < 0)
{
// Imaginary roots only - face not crossed
return GREAT;
}
else
{
scalar q = -0.5*(b + sign(b)*Foam::sqrt(discriminant));
if (mag(q) < VSMALL)
{
// If q is zero, then l1 = q/a is the required
// value of lambda, and is zero.
return 0.0;
}
scalar l1 = q/a;
scalar l2 = c/q;
// There will be two roots, a big one and a little
// one, choose the little one.
if (mag(l1) < mag(l2))
{
return l1;
}
else
{
return l2;
}
}
}
{
// When a is zero, solve the first order polynomial
lambdaNumerator = -c;
lambdaDenominator = b;
}
}
else
{
// When n = n0 is zero, there is no plane rotation, solve the
// first order polynomial
lambdaNumerator = -(dS & n0);
lambdaDenominator = ((dP - dB) & n0);
}
if (mag(lambdaDenominator) < tol)
{
if (mag(lambdaNumerator) < tol)
{
// Track starts on the face, and is potentially
// parallel to it. +-tol)/+-tol is not a good
// comparison, return 0.0, in anticipation of tet
// centre correction.
return 0.0;
}
else
{
if (mag((to - from)) < tol/mag(n))
{
// Zero length track, not along the face, face
// cannot be crossed.
return GREAT;
}
else
{
// Trajectory is non-zero and parallel to face
lambdaDenominator = sign(lambdaDenominator)*SMALL;
}
}
}
return lambdaNumerator/lambdaDenominator;
}
inline void Foam::particle::tetNeighbour(label triI)
{
const labelList& pOwner = mesh_.faceOwner();
const faceList& pFaces = mesh_.faces();
bool own = (pOwner[tetFacei_] == celli_);
const Foam::face& f = pFaces[tetFacei_];
label tetBasePtI = mesh_.tetBasePtIs()[tetFacei_];
if (tetBasePtI == -1)
{
FatalErrorInFunction
<< "No base point for face " << tetFacei_ << ", " << f
<< ", produces a valid tet decomposition."
<< abort(FatalError);
}
label facePtI = (tetPti_ + tetBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
switch (triI)
{
case 0:
{
// Crossing this triangle changes tet to that in the
// neighbour cell over tetFacei
// Modification of celli_ will happen by other indexing,
// tetFacei_ and tetPti don't change.
break;
}
case 1:
{
crossEdgeConnectedFace
(
celli_,
tetFacei_,
tetPti_,
Foam::edge(f[facePtI], f[otherFacePtI])
);
break;
}
case 2:
{
if (own)
{
if (tetPti_ < f.size() - 2)
{
tetPti_ = f.fcIndex(tetPti_);
}
else
{
crossEdgeConnectedFace
(
celli_,
tetFacei_,
tetPti_,
Foam::edge(f[tetBasePtI], f[otherFacePtI])
);
}
}
else
{
if (tetPti_ > 1)
{
tetPti_ = f.rcIndex(tetPti_);
}
else
{
crossEdgeConnectedFace
(
celli_,
tetFacei_,
tetPti_,
Foam::edge(f[tetBasePtI], f[facePtI])
);
}
}
break;
}
case 3:
{
if (own)
{
if (tetPti_ > 1)
{
tetPti_ = f.rcIndex(tetPti_);
}
else
{
crossEdgeConnectedFace
(
celli_,
tetFacei_,
tetPti_,
Foam::edge(f[tetBasePtI], f[facePtI])
);
}
}
else
{
if (tetPti_ < f.size() - 2)
{
tetPti_ = f.fcIndex(tetPti_);
}
else
{
crossEdgeConnectedFace
(
celli_,
tetFacei_,
tetPti_,
Foam::edge(f[tetBasePtI], f[otherFacePtI])
);
}
}
break;
}
default:
{
FatalErrorInFunction
<< "Tet tri face index error, can only be 0..3, supplied "
<< triI << abort(FatalError);
break;
}
}
}
inline void Foam::particle::crossEdgeConnectedFace
(
const label& celli,
label& tetFacei,
label& tetPti,
const edge& e
)
{
const faceList& pFaces = mesh_.faces();
const cellList& pCells = mesh_.cells();
const Foam::face& f = pFaces[tetFacei];
const Foam::cell& thisCell = pCells[celli];
forAll(thisCell, cFI)
{
// Loop over all other faces of this cell and
// find the one that shares this edge
label fI = thisCell[cFI];
if (tetFacei == fI)
{
continue;
}
const Foam::face& otherFace = pFaces[fI];
label edDir = otherFace.edgeDirection(e);
if (edDir == 0)
{
continue;
}
else if (f == pFaces[fI])
{
// This is a necessary condition if using duplicate baffles
// (so coincident faces). We need to make sure we don't cross into
// the face with the same vertices since we might enter a tracking
// loop where it never exits. This test should be cheap
// for most meshes so can be left in for 'normal' meshes.
continue;
}
else
{
//Found edge on other face
tetFacei = fI;
label eIndex = -1;
if (edDir == 1)
{
// Edge is in the forward circulation of this face, so
// work with the start point of the edge
eIndex = findIndex(otherFace, e.start());
}
else
{
// edDir == -1, so the edge is in the reverse
// circulation of this face, so work with the end
// point of the edge
eIndex = findIndex(otherFace, e.end());
}
label tetBasePtI = mesh_.tetBasePtIs()[fI];
if (tetBasePtI == -1)
{
FatalErrorInFunction
<< "No base point for face " << fI << ", " << f
<< ", produces a decomposition that has a minimum "
<< "volume greater than tolerance."
<< abort(FatalError);
}
// Find eIndex relative to the base point on new face
eIndex -= tetBasePtI;
if (neg(eIndex))
{
eIndex = (eIndex + otherFace.size()) % otherFace.size();
}
if (eIndex == 0)
{
// The point is the base point, so this is first tet
// in the face circulation
tetPti = 1;
}
else if (eIndex == otherFace.size() - 1)
{
// The point is the last before the base point, so
// this is the last tet in the face circulation
tetPti = otherFace.size() - 2;
}
else
{
tetPti = eIndex;
}
break;
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::particle::getNewParticleID() const inline Foam::label Foam::particle::getNewParticleID() const
@ -583,15 +48,9 @@ inline const Foam::polyMesh& Foam::particle::mesh() const
} }
inline const Foam::vector& Foam::particle::position() const inline const Foam::barycentric& Foam::particle::coordinates() const
{ {
return position_; return coordinates_;
}
inline Foam::vector& Foam::particle::position()
{
return position_;
} }
@ -601,36 +60,18 @@ inline Foam::label Foam::particle::cell() const
} }
inline Foam::label& Foam::particle::cell()
{
return celli_;
}
inline Foam::label Foam::particle::tetFace() const inline Foam::label Foam::particle::tetFace() const
{ {
return tetFacei_; return tetFacei_;
} }
inline Foam::label& Foam::particle::tetFace()
{
return tetFacei_;
}
inline Foam::label Foam::particle::tetPt() const inline Foam::label Foam::particle::tetPt() const
{ {
return tetPti_; return tetPti_;
} }
inline Foam::label& Foam::particle::tetPt()
{
return tetPti_;
}
inline Foam::tetIndices Foam::particle::currentTetIndices() const inline Foam::tetIndices Foam::particle::currentTetIndices() const
{ {
return tetIndices(celli_, tetFacei_, tetPti_, mesh_); return tetIndices(celli_, tetFacei_, tetPti_, mesh_);
@ -661,161 +102,37 @@ inline Foam::label Foam::particle::face() const
} }
inline Foam::label& Foam::particle::face() inline bool Foam::particle::onFace() const
{ {
return facei_; return facei_ >= 0;
} }
inline void Foam::particle::initCellFacePt() inline bool Foam::particle::onInternalFace() const
{ {
if (celli_ == -1) return onFace() && mesh_.isInternalFace(facei_);
{ }
mesh_.findCellFacePt
(
position_,
celli_,
tetFacei_,
tetPti_
);
if (debug && celli_ == -1)
{ inline bool Foam::particle::onBoundaryFace() const
WarningInFunction {
<< "cell, tetFace and tetPt search failure for position " return onFace() && !mesh_.isInternalFace(facei_);
<< position_ }
<< endl;
}
inline Foam::vector Foam::particle::position() const
{
if (mesh_.moving())
{
return movingTetTransform(0)[0] & coordinates_;
} }
else else
{ {
mesh_.findTetFacePt(celli_, position_, tetFacei_, tetPti_); return tetTransform() & coordinates_;
if (tetFacei_ == -1 || tetPti_ == -1)
{
label oldCelli = celli_;
mesh_.findCellFacePt
(
position_,
celli_,
tetFacei_,
tetPti_
);
if (celli_ == -1 || tetFacei_ == -1 || tetPti_ == -1)
{
// The particle has entered this function with a cell number,
// but hasn't been able to find a cell to occupy.
if (!mesh_.pointInCellBB(position_, oldCelli, 0.1))
{
// If the position is not inside the (slightly extended)
// bound-box of the cell that it thought it should be in,
// then this is considered an error.
if (debug)
{
WarningInFunction
<< "position " << position_
<< " not in specified cell " << oldCelli
<< " with centre " << mesh_.cellCentres()[oldCelli]
<< endl;
}
celli_ = -1;
tetFacei_ = -1;
tetPti_ = -1;
return;
}
// The position is in the (slightly extended) bound-box of the
// cell. This situation may arise because the face
// decomposition of the cell is not the same as when the
// particle acquired the cell index. For example, it has been
// read into a mesh that has made a different face base-point
// decision for a boundary face and now this particle is in a
// position that is not in the mesh. Gradually move the
// particle towards the centre of the cell that it thought that
// it was in.
celli_ = oldCelli;
point newPosition = position_;
const point& cC = mesh_.cellCentres()[celli_];
label trap(1.0/trackingCorrectionTol + 1);
label iterNo = 0;
do
{
newPosition += trackingCorrectionTol*(cC - position_);
mesh_.findTetFacePt
(
celli_,
newPosition,
tetFacei_,
tetPti_
);
iterNo++;
} while (tetFacei_ < 0 && iterNo <= trap);
if (tetFacei_ == -1)
{
FatalErrorInFunction
<< "cell, tetFace and tetPt search failure at position "
<< position_ << abort(FatalError);
}
if (debug)
{
WarningInFunction
<< "Particle moved from " << position_
<< " to " << newPosition
<< " in cell " << celli_
<< " tetFace " << tetFacei_
<< " tetPt " << tetPti_ << nl
<< " (A fraction of "
<< 1.0 - mag(cC - newPosition)/mag(cC - position_)
<< " of the distance to the cell centre)"
<< " because a decomposition tetFace and tetPt "
<< "could not be found."
<< endl;
}
position_ = newPosition;
}
if (debug && celli_ != oldCelli)
{
WarningInFunction
<< "Particle at position " << position_
<< " searched for a cell, tetFace and tetPt." << nl
<< " Found"
<< " cell " << celli_
<< " tetFace " << tetFacei_
<< " tetPt " << tetPti_ << nl
<< " This is a different cell to that which was supplied"
<< " (" << oldCelli << ")." << nl
<< endl;
}
}
} }
} }
inline bool Foam::particle::onBoundary() const
{
return facei_ != -1 && facei_ >= mesh_.nInternalFaces();
}
inline Foam::scalar& Foam::particle::stepFraction() inline Foam::scalar& Foam::particle::stepFraction()
{ {
return stepFraction_; return stepFraction_;
@ -852,29 +169,9 @@ inline Foam::label& Foam::particle::origId()
} }
inline bool Foam::particle::softImpact() const
{
return false;
}
inline Foam::scalar Foam::particle::currentTime() const inline Foam::scalar Foam::particle::currentTime() const
{ {
return return mesh_.time().value() + stepFraction_*mesh_.time().deltaTValue();
mesh_.time().value()
+ stepFraction_*mesh_.time().deltaTValue();
}
inline bool Foam::particle::internalFace(const label facei) const
{
return mesh_.isInternalFace(facei);
}
bool Foam::particle::boundaryFace(const label facei) const
{
return !internalFace(facei);
} }
@ -894,10 +191,4 @@ inline Foam::label Foam::particle::patchFace
} }
inline Foam::label Foam::particle::faceInterpolation() const
{
return facei_;
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -32,12 +32,12 @@ Foam::string Foam::particle::propertyList_ = Foam::particle::propertyList();
const std::size_t Foam::particle::sizeofPosition_ const std::size_t Foam::particle::sizeofPosition_
( (
offsetof(particle, facei_) - offsetof(particle, position_) offsetof(particle, facei_) - offsetof(particle, coordinates_)
); );
const std::size_t Foam::particle::sizeofFields_ const std::size_t Foam::particle::sizeofFields_
( (
sizeof(particle) - offsetof(particle, position_) sizeof(particle) - offsetof(particle, coordinates_)
); );
@ -46,38 +46,33 @@ const std::size_t Foam::particle::sizeofFields_
Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields) Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
: :
mesh_(mesh), mesh_(mesh),
position_(), coordinates_(),
celli_(-1), celli_(-1),
facei_(-1),
stepFraction_(0.0),
tetFacei_(-1), tetFacei_(-1),
tetPti_(-1), tetPti_(-1),
facei_(-1),
stepFraction_(0.0),
origProc_(Pstream::myProcNo()), origProc_(Pstream::myProcNo()),
origId_(-1) origId_(-1)
{ {
if (is.format() == IOstream::ASCII) if (is.format() == IOstream::ASCII)
{ {
is >> position_ >> celli_; is >> coordinates_ >> celli_ >> tetFacei_ >> tetPti_;
if (readFields) if (readFields)
{ {
is >> facei_ is >> facei_ >> stepFraction_ >> origProc_ >> origId_;
>> stepFraction_
>> tetFacei_
>> tetPti_
>> origProc_
>> origId_;
} }
} }
else else
{ {
if (readFields) if (readFields)
{ {
is.read(reinterpret_cast<char*>(&position_), sizeofFields_); is.read(reinterpret_cast<char*>(&coordinates_), sizeofFields_);
} }
else else
{ {
is.read(reinterpret_cast<char*>(&position_), sizeofPosition_); is.read(reinterpret_cast<char*>(&coordinates_), sizeofPosition_);
} }
} }
@ -90,11 +85,14 @@ void Foam::particle::writePosition(Ostream& os) const
{ {
if (os.format() == IOstream::ASCII) if (os.format() == IOstream::ASCII)
{ {
os << position_ << token::SPACE << celli_; os << coordinates_
<< token::SPACE << celli_
<< token::SPACE << tetFacei_
<< token::SPACE << tetPti_;
} }
else else
{ {
os.write(reinterpret_cast<const char*>(&position_), sizeofPosition_); os.write(reinterpret_cast<const char*>(&coordinates_), sizeofPosition_);
} }
// Check state of Ostream // Check state of Ostream
@ -106,12 +104,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
{ {
if (os.format() == IOstream::ASCII) if (os.format() == IOstream::ASCII)
{ {
os << p.position_ os << p.coordinates_
<< token::SPACE << p.celli_ << token::SPACE << p.celli_
<< token::SPACE << p.facei_
<< token::SPACE << p.stepFraction_
<< token::SPACE << p.tetFacei_ << token::SPACE << p.tetFacei_
<< token::SPACE << p.tetPti_ << token::SPACE << p.tetPti_
<< token::SPACE << p.facei_
<< token::SPACE << p.stepFraction_
<< token::SPACE << p.origProc_ << token::SPACE << p.origProc_
<< token::SPACE << p.origId_; << token::SPACE << p.origId_;
} }
@ -119,7 +117,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
{ {
os.write os.write
( (
reinterpret_cast<const char*>(&p.position_), reinterpret_cast<const char*>(&p.coordinates_),
particle::sizeofFields_ particle::sizeofFields_
); );
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -44,7 +44,7 @@ void Foam::particle::prepareForParallelTransfer
) )
{ {
// Convert the face index to be local to the processor patch // Convert the face index to be local to the processor patch
facei_ = patchFace(patchi, facei_); facei_ = mesh_.boundaryMesh()[patchi].whichFace(facei_);
} }
@ -58,12 +58,6 @@ void Foam::particle::correctAfterParallelTransfer
const coupledPolyPatch& ppp = const coupledPolyPatch& ppp =
refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]); refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
celli_ = ppp.faceCells()[facei_];
// Have patch transform the position
ppp.transformPosition(position_, facei_);
// Transform the properties
if (!ppp.parallel()) if (!ppp.parallel())
{ {
const tensor& T = const tensor& T =
@ -85,41 +79,22 @@ void Foam::particle::correctAfterParallelTransfer
transformProperties(-s); transformProperties(-s);
} }
tetFacei_ = facei_ + ppp.start(); // Set the topology
celli_ = ppp.faceCells()[facei_];
// Faces either side of a coupled patch have matched base indices, facei_ += ppp.start();
// tetPti is specified relative to the base point, already and tetFacei_ = facei_;
// opposite circulation directions by design, so if the vertices // Faces either side of a coupled patch are numbered in opposite directions
// are: // as their normals both point away from their connected cells. The tet
// source: // point therefore counts in the opposite direction from the base point.
// face (a b c d e f)
// fPtI 0 1 2 3 4 5
// +
// destination:
// face (a f e d c b)
// fPtI 0 1 2 3 4 5
// +
// where a is the base point of the face are matching , and we
// have fPtI = 1 on the source processor face, i.e. vertex b, then
// this because of the face circulation direction change, vertex c
// is the characterising point on the destination processor face,
// giving the destination fPtI as:
// fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
// This relationship can be verified for other points and sizes of
// face.
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_; tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
// Reset the face index for the next tracking operation // Reflect to account for the change of triangle orientation in the new cell
if (stepFraction_ > (1.0 - SMALL)) reflect();
{
stepFraction_ = 1.0; // Note that the position does not need transforming explicitly. The face-
facei_ = -1; // triangle on the receive patch is the transformation of the one on the
} // send patch, so whilst the barycentric coordinates remain the same, the
else // change of triangle implicitly transforms the position.
{
facei_ += ppp.start();
}
} }
@ -158,18 +133,20 @@ void Foam::particle::readFields(CloudType& c)
template<class CloudType> template<class CloudType>
void Foam::particle::writeFields(const CloudType& c) void Foam::particle::writeFields(const CloudType& c)
{ {
// Write the cloud position file
IOPosition<CloudType> ioP(c); IOPosition<CloudType> ioP(c);
ioP.write(); ioP.write();
label np = c.size(); label np = c.size();
IOField<label> origProc IOField<label> origProc
( (
c.fieldIOobject("origProcId", IOobject::NO_READ), c.fieldIOobject("origProcId", IOobject::NO_READ),
np np
); );
IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np); IOField<label> origId
(
c.fieldIOobject("origId", IOobject::NO_READ),
np
);
label i = 0; label i = 0;
forAllConstIter(typename CloudType, c, iter) forAllConstIter(typename CloudType, c, iter)
@ -185,387 +162,55 @@ void Foam::particle::writeFields(const CloudType& c)
template<class TrackData> template<class TrackData>
Foam::label Foam::particle::track(const vector& endPosition, TrackData& td) void Foam::particle::trackToFace
{
facei_ = -1;
// Tracks to endPosition or stop on boundary
while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
{
stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
}
return facei_;
}
template<class TrackData>
Foam::scalar Foam::particle::trackToFace
( (
const vector& endPosition, const vector& displacement,
const scalar fraction,
TrackData& td TrackData& td
) )
{ {
typedef typename TrackData::cloudType cloudType; // Track
typedef typename cloudType::particleType particleType; trackToFace(displacement, fraction);
cloudType& cloud = td.cloud(); // If the track is complete, return
if (!onFace())
const faceList& pFaces = mesh_.faces();
const pointField& pPts = mesh_.points();
const vectorField& pC = mesh_.cellCentres();
facei_ = -1;
// Pout<< "Particle " << origId_ << " " << origProc_
// << " Tracking from " << position_
// << " to " << endPosition
// << endl;
// Pout<< "stepFraction " << stepFraction_ << nl
// << "celli " << celli_ << nl
// << "tetFacei " << tetFacei_ << nl
// << "tetPti " << tetPti_
// << endl;
scalar trackFraction = 0.0;
// Minimum tetrahedron decomposition of each cell of the mesh into
// using the cell centre, base point on face, and further two
// points on the face. For each face of n points, there are n - 2
// tets generated.
// The points for each tet are organised to match those used in the
// tetrahedron class, supplying them in the order:
// Cc, basePt, pA, pB
// where:
// + Cc is the cell centre;
// + basePt is the base point on the face;
// + pA and pB are the remaining points on the face, such that
// the circulation, {basePt, pA, pB} produces a positive
// normal by the right-hand rule. pA and pB are chosen from
// tetPti_ do accomplish this depending if the cell owns the
// face, tetPti_ is the vertex that characterises the tet, and
// is the first vertex on the tet when circulating around the
// face. Therefore, the same tetPti represents the same face
// triangle for both the owner and neighbour cell.
//
// Each tet has its four triangles represented in the same order:
// 0) tri joining a tet to the tet across the face in next cell.
// This is the triangle opposite Cc.
// 1) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the edge of the current face
// that doesn't contain basePt. This is the triangle opposite
// basePt.
// 2) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the tet-edge (basePt - pB).
// This may be on the same face, or a different one. This is
// the triangle opposite basePt. This is the triangle opposite
// pA.
// 4) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the tet-edge (basePt - pA).
// This may be on the same face, or a different one. This is
// the triangle opposite basePt. This is the triangle opposite
// pA.
// Which tri (0..3) of the tet has been crossed
label triI = -1;
// Determine which face was actually crossed. lambdaMin < SMALL
// is considered a trigger for a tracking correction towards the
// current tet centre.
scalar lambdaMin = VGREAT;
DynamicList<label>& tris = cloud.labels();
// Tet indices that will be set by hitWallFaces if a wall face is
// to be hit, or are set when any wall tri of a tet is hit.
// Carries the description of the tet on which the cell face has
// been hit. For the case of being set in hitWallFaces, this may
// be a different tet to the one that the particle occupies.
tetIndices faceHitTetIs;
// What tolerance is appropriate the minimum lambda numerator and
// denominator for tracking in this cell.
scalar lambdaDistanceTolerance =
lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[celli_];
do
{ {
if (triI != -1) return;
{ }
// Change tet ownership because a tri face has been crossed
tetNeighbour(triI);
}
const Foam::face& f = pFaces[tetFacei_];
bool own = (mesh_.faceOwner()[tetFacei_] == celli_);
label tetBasePtI = mesh_.tetBasePtIs()[tetFacei_];
label basePtI = f[tetBasePtI];
label facePtI = (tetPti_ + tetBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label fPtAI = -1;
label fPtBI = -1;
if (own)
{
fPtAI = facePtI;
fPtBI = otherFacePtI;
}
else
{
fPtAI = otherFacePtI;
fPtBI = facePtI;
}
tetPointRef tet
(
pC[celli_],
pPts[basePtI],
pPts[f[fPtAI]],
pPts[f[fPtBI]]
);
if (lambdaMin < SMALL)
{
// Apply tracking correction towards tet centre
if (debug)
{
Pout<< "tracking rescue using tetCentre from " << position();
}
position_ += trackingCorrectionTol*(tet.centre() - position_);
if (debug)
{
Pout<< " to " << position() << " due to "
<< (tet.centre() - position_) << endl;
}
cloud.trackingRescue();
return trackFraction;
}
if (triI != -1 && mesh_.moving())
{
// Mesh motion requires stepFraction to be correct for
// each tracking portion, so trackToFace must return after
// every lambda calculation.
return trackFraction;
}
FixedList<vector, 4> tetAreas;
tetAreas[0] = tet.Sa();
tetAreas[1] = tet.Sb();
tetAreas[2] = tet.Sc();
tetAreas[3] = tet.Sd();
FixedList<label, 4> tetPlaneBasePtIs;
tetPlaneBasePtIs[0] = basePtI;
tetPlaneBasePtIs[1] = f[fPtAI];
tetPlaneBasePtIs[2] = basePtI;
tetPlaneBasePtIs[3] = basePtI;
findTris
(
endPosition,
tris,
tet,
tetAreas,
tetPlaneBasePtIs,
lambdaDistanceTolerance
);
// Reset variables for new track
triI = -1;
lambdaMin = VGREAT;
// Pout<< "tris " << tris << endl;
// Sets a value for lambdaMin and facei_ if a wall face is hit
// by the track.
hitWallFaces
(
cloud,
position_,
endPosition,
lambdaMin,
faceHitTetIs
);
// Did not hit any tet tri faces, and no wall face has been
// found to hit.
if (tris.empty() && facei_ < 0)
{
position_ = endPosition;
return 1.0;
}
else
{
// Loop over all found tris and see if any of them find a
// lambda value smaller than that found for a wall face.
forAll(tris, i)
{
label tI = tris[i];
scalar lam = tetLambda
(
position_,
endPosition,
tI,
tetAreas[tI],
tetPlaneBasePtIs[tI],
celli_,
tetFacei_,
tetPti_,
lambdaDistanceTolerance
);
if (lam < lambdaMin)
{
lambdaMin = lam;
triI = tI;
}
}
}
if (triI == 0)
{
// This must be a cell face crossing
facei_ = tetFacei_;
// Set the faceHitTetIs to those for the current tet in case a
// wall interaction is required with the cell face
faceHitTetIs = tetIndices
(
celli_,
tetFacei_,
tetBasePtI,
fPtAI,
fPtBI,
tetPti_
);
}
else if (triI > 0)
{
// A tri was found to be crossed before a wall face was hit (if any)
facei_ = -1;
}
// Pout<< "track loop " << position_ << " " << endPosition << nl
// << " " << celli_
// << " " << facei_
// << " " << tetFacei_
// << " " << tetPti_
// << " " << triI
// << " " << lambdaMin
// << " " << trackFraction
// << endl;
// Pout<< "# Tracking loop tet "
// << origId_ << " " << origProc_<< nl
// << "# face: " << tetFacei_ << nl
// << "# tetPti: " << tetPti_ << nl
// << "# tetBasePtI: " << mesh_.tetBasePtIs()[tetFacei_] << nl
// << "# tet.mag(): " << tet.mag() << nl
// << "# tet.quality(): " << tet.quality()
// << endl;
// meshTools::writeOBJ(Pout, tet.a());
// meshTools::writeOBJ(Pout, tet.b());
// meshTools::writeOBJ(Pout, tet.c());
// meshTools::writeOBJ(Pout, tet.d());
// Pout<< "f 1 3 2" << nl
// << "f 2 3 4" << nl
// << "f 1 4 3" << nl
// << "f 1 2 4" << endl;
// The particle can be 'outside' the tet. This will yield a
// lambda larger than 1, or smaller than 0. For values < 0,
// the particle travels away from the tet and we don't move
// the particle, only change tet/cell. For values larger than
// 1, we move the particle to endPosition before the tet/cell
// change.
if (lambdaMin > SMALL)
{
if (lambdaMin <= 1.0)
{
trackFraction += lambdaMin*(1 - trackFraction);
position_ += lambdaMin*(endPosition - position_);
}
else
{
position_ = endPosition;
return 1.0;
}
}
else
{
// Set lambdaMin to zero to force a towards-tet-centre
// correction.
lambdaMin = 0.0;
}
} while (facei_ < 0);
// Hit face/patch processing
typedef typename TrackData::cloudType::particleType particleType;
particleType& p = static_cast<particleType&>(*this); particleType& p = static_cast<particleType&>(*this);
p.hitFace(td); p.hitFace(td);
if (onInternalFace())
if (internalFace(facei_))
{ {
// Change tet ownership because a tri face has been crossed, changeCell();
// in general this is:
// tetNeighbour(triI);
// but triI must be 0;
// No modifications are required for triI = 0, no call required to
// tetNeighbour(0);
if (celli_ == mesh_.faceOwner()[facei_])
{
celli_ = mesh_.faceNeighbour()[facei_];
}
else if (celli_ == mesh_.faceNeighbour()[facei_])
{
celli_ = mesh_.faceOwner()[facei_];
}
else
{
FatalErrorInFunction
<< "addressing failure" << abort(FatalError);
}
} }
else else
{ {
label origFacei = facei_; label origFacei = facei_;
label patchi = patch(facei_); label patchi = mesh_.boundaryMesh().whichPatch(facei_);
// No action taken for tetPti_ for tetFacei_ here, handled by // No action is taken for tetPti_ for tetFacei_ here. These are handled
// patch interaction call or later during processor transfer. // by the patch interaction call or later during processor transfer.
const tetIndices faceHitTetIs =
polyMeshTetDecomposition::triangleTetIndices
(
mesh_,
tetFacei_,
celli_,
tetPti_
);
if if
( (
!p.hitPatch !p.hitPatch
( (
mesh_.boundaryMesh()[patchi], mesh_.boundaryMesh()[patchi],
td, td,
patchi, patchi,
trackFraction, stepFraction(),
faceHitTetIs faceHitTetIs
) )
) )
@ -573,7 +218,7 @@ Foam::scalar Foam::particle::trackToFace
// Did patch interaction model switch patches? // Did patch interaction model switch patches?
if (facei_ != origFacei) if (facei_ != origFacei)
{ {
patchi = patch(facei_); patchi = mesh_.boundaryMesh().whichPatch(facei_);
} }
const polyPatch& patch = mesh_.boundaryMesh()[patchi]; const polyPatch& patch = mesh_.boundaryMesh()[patchi];
@ -612,7 +257,7 @@ Foam::scalar Foam::particle::trackToFace
( (
static_cast<const cyclicAMIPolyPatch&>(patch), static_cast<const cyclicAMIPolyPatch&>(patch),
td, td,
endPosition - position_ displacement
); );
} }
else if (isA<processorPolyPatch>(patch)) else if (isA<processorPolyPatch>(patch))
@ -635,290 +280,6 @@ Foam::scalar Foam::particle::trackToFace
} }
} }
} }
if (lambdaMin < SMALL)
{
// Apply tracking correction towards tet centre.
// Generate current tet to find centre to apply correction.
tetPointRef tet = currentTet();
if (debug)
{
Pout<< "tracking rescue for lambdaMin:" << lambdaMin
<< "from " << position();
}
position_ += trackingCorrectionTol*(tet.centre() - position_);
if
(
cloud.hasWallImpactDistance()
&& !internalFace(faceHitTetIs.face())
&& cloud.cellHasWallFaces()[faceHitTetIs.cell()]
)
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
label fI = faceHitTetIs.face();
label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
if (isA<wallPolyPatch>(patches[patchi]))
{
// In the case of collision with a wall where there is
// a non-zero wallImpactDistance, it is possible for
// there to be a tracking correction required to bring
// the particle into the domain, but the position of
// the particle is further from the wall than the tet
// centre, in which case the normal correction can be
// counter-productive, i.e. pushes the particle
// further out of the domain. In this case it is the
// position that hit the wall that is in need of a
// rescue correction.
triPointRef wallTri = faceHitTetIs.faceTri(mesh_);
tetPointRef wallTet = faceHitTetIs.tet(mesh_);
vector nHat = wallTri.normal();
nHat /= mag(nHat);
const scalar r = p.wallImpactDistance(nHat);
// Removing (approximately) the wallTri normal
// component of the existing correction, to avoid the
// situation where the existing correction in the wall
// normal direction is larger towards the wall than
// the new correction is away from it.
position_ +=
trackingCorrectionTol
*(
(wallTet.centre() - (position_ + r*nHat))
- (nHat & (tet.centre() - position_))*nHat
);
}
}
if (debug)
{
Pout<< " to " << position() << endl;
}
cloud.trackingRescue();
}
return trackFraction;
}
template<class CloudType>
void Foam::particle::hitWallFaces
(
const CloudType& cloud,
const vector& from,
const vector& to,
scalar& lambdaMin,
tetIndices& closestTetIs
)
{
typedef typename CloudType::particleType particleType;
if (!(cloud.hasWallImpactDistance() && cloud.cellHasWallFaces()[celli_]))
{
return;
}
particleType& p = static_cast<particleType&>(*this);
const faceList& pFaces = mesh_.faces();
const Foam::cell& thisCell = mesh_.cells()[celli_];
scalar lambdaDistanceTolerance =
lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[celli_];
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(thisCell, cFI)
{
label fI = thisCell[cFI];
if (internalFace(fI))
{
continue;
}
label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
if (isA<wallPolyPatch>(patches[patchi]))
{
// Get the decomposition of this wall face
const List<tetIndices> faceTetIs =
polyMeshTetDecomposition::faceTetIndices(mesh_, fI, celli_);
const Foam::face& f = pFaces[fI];
forAll(faceTetIs, tI)
{
const tetIndices& tetIs = faceTetIs[tI];
triPointRef tri = tetIs.faceTri(mesh_);
vector n = tri.normal();
vector nHat = n/mag(n);
// Radius of particle with respect to this wall face
// triangle. Assuming that the wallImpactDistance
// does not change as the particle or the mesh moves
// forward in time.
scalar r = p.wallImpactDistance(nHat);
vector toPlusRNHat = to + r*nHat;
// triI = 0 because it is the cell face tri of the tet
// we are concerned with.
scalar tetClambda = tetLambda
(
tetIs.tet(mesh_).centre(),
toPlusRNHat,
0,
n,
f[tetIs.faceBasePt()],
celli_,
fI,
tetIs.tetPt(),
lambdaDistanceTolerance
);
if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
{
// toPlusRNHat is not on the outside of the plane of
// the wall face tri, the tri cannot be hit.
continue;
}
// Check if the actual trajectory of the near-tri
// points intersects the triangle.
vector fromPlusRNHat = from + r*nHat;
// triI = 0 because it is the cell face tri of the tet
// we are concerned with.
scalar lambda = tetLambda
(
fromPlusRNHat,
toPlusRNHat,
0,
n,
f[tetIs.faceBasePt()],
celli_,
fI,
tetIs.tetPt(),
lambdaDistanceTolerance
);
pointHit hitInfo(Zero);
if (mesh_.moving())
{
// For a moving mesh, the position of wall
// triangle needs to be moved in time to be
// consistent with the moment defined by the
// current value of stepFraction and the value of
// lambda just calculated.
// Total fraction thought the timestep of the
// motion, including stepFraction before the
// current tracking step and the current
// lambda
// i.e.
// let s = stepFraction, l = lambda
// Motion of x in time:
// |-----------------|---------|---------|
// x00 x0 xi x
//
// where xi is the correct value of x at the required
// tracking instant.
//
// x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
//
// i.e. the motion covered by previous tracking portions
// within this timestep, and
//
// xi = x0 + l*(x - x0)
// = l*x + (1 - l)*x0
// = l*x + (1 - l)*(s*x + (1 - s)*x00)
// = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
//
// let m = (s + l - s*l)
//
// xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
//
// In the same form as before.
// Clip lambda to 0.0-1.0 to ensure that sensible
// positions are used for triangle intersections.
scalar lam = max(0.0, min(1.0, lambda));
scalar m = stepFraction_ + lam - (stepFraction_*lam);
triPointRef tri00 = tetIs.oldFaceTri(mesh_);
// Use SMALL positive tolerance to make the triangle
// slightly "fat" to improve robustness. Intersection
// is calculated as the ray (from + r*nHat) -> (to +
// r*nHat).
point tPtA = tri00.a() + m*(tri.a() - tri00.a());
point tPtB = tri00.b() + m*(tri.b() - tri00.b());
point tPtC = tri00.c() + m*(tri.c() - tri00.c());
triPointRef t(tPtA, tPtB, tPtC);
// The point fromPlusRNHat + m*(to - from) is on the
// plane of the triangle. Determine the
// intersection with this triangle by testing if
// this point is inside or outside of the triangle.
hitInfo = t.intersection
(
fromPlusRNHat + m*(to - from),
t.normal(),
intersection::FULL_RAY,
SMALL
);
}
else
{
// Use SMALL positive tolerance to make the triangle
// slightly "fat" to improve robustness. Intersection
// is calculated as the ray (from + r*nHat) -> (to +
// r*nHat).
hitInfo = tri.intersection
(
fromPlusRNHat,
(to - from),
intersection::FULL_RAY,
SMALL
);
}
if (hitInfo.hit())
{
if (lambda < lambdaMin)
{
lambdaMin = lambda;
facei_ = fI;
closestTetIs = tetIs;
}
}
}
}
}
} }
@ -994,22 +355,17 @@ void Foam::particle::hitCyclicPatch
TrackData& td TrackData& td
) )
{ {
facei_ = cpp.transformGlobalFace(facei_); const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
const label receiveFacei = receiveCpp.whichFace(facei_);
// Set the topology
facei_ = tetFacei_ = cpp.transformGlobalFace(facei_);
celli_ = mesh_.faceOwner()[facei_]; celli_ = mesh_.faceOwner()[facei_];
// See note in correctAfterParallelTransfer for tetPti addressing ...
tetFacei_ = facei_;
// See note in correctAfterParallelTransfer for tetPti_ addressing.
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_; tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
const cyclicPolyPatch& receiveCpp = cpp.neighbPatch(); // Reflect to account for the change of triangle orientation in the new cell
label patchFacei = receiveCpp.whichFace(facei_); reflect();
// Now the particle is on the receiving side
// Have patch transform the position
receiveCpp.transformPosition(position_, patchFacei);
// Transform the properties // Transform the properties
if (!receiveCpp.parallel()) if (!receiveCpp.parallel())
@ -1018,7 +374,7 @@ void Foam::particle::hitCyclicPatch
( (
receiveCpp.forwardT().size() == 1 receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0] ? receiveCpp.forwardT()[0]
: receiveCpp.forwardT()[patchFacei] : receiveCpp.forwardT()[receiveFacei]
); );
transformProperties(T); transformProperties(T);
} }
@ -1028,7 +384,7 @@ void Foam::particle::hitCyclicPatch
( (
(receiveCpp.separation().size() == 1) (receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0] ? receiveCpp.separation()[0]
: receiveCpp.separation()[patchFacei] : receiveCpp.separation()[receiveFacei]
); );
transformProperties(-s); transformProperties(-s);
} }
@ -1043,38 +399,41 @@ void Foam::particle::hitCyclicAMIPatch
const vector& direction const vector& direction
) )
{ {
vector pos = position();
const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch(); const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
const label sendFacei = cpp.whichFace(facei_);
const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
// Patch face index on sending side if (receiveFacei < 0)
label patchFacei = facei_ - cpp.start();
// Patch face index on receiving side - also updates position
patchFacei = cpp.pointFace(patchFacei, direction, position_);
if (patchFacei < 0)
{ {
// If the patch face of the particle is not known assume that // If the patch face of the particle is not known assume that the
// the particle is lost and to be deleted // particle is lost and mark it to be deleted.
td.keepParticle = false; td.keepParticle = false;
WarningInFunction WarningInFunction
<< "Particle lost across " << cyclicAMIPolyPatch::typeName << "Particle lost across " << cyclicAMIPolyPatch::typeName
<< " patches " << cpp.name() << " and " << receiveCpp.name() << " patches " << cpp.name() << " and " << receiveCpp.name()
<< " at position " << position_ << " at position " << pos << endl;
<< endl;
} }
// Convert face index into global numbering // Set the topology
facei_ = patchFacei + receiveCpp.start(); facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
celli_ = mesh_.faceOwner()[facei_]; // Locate the particle on the recieving side
locate
(
pos,
&direction,
mesh_.faceOwner()[facei_],
false,
"Particle crossed between " + cyclicAMIPolyPatch::typeName +
" patches " + cpp.name() + " and " + receiveCpp.name() +
" to a location outside of the mesh."
);
tetFacei_ = facei_; // The particle must remain associated with a face for the tracking to
// register as incomplete
// See note in correctAfterParallelTransfer for tetPti_ addressing. facei_ = tetFacei_;
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
// Now the particle is on the receiving side
// Transform the properties // Transform the properties
if (!receiveCpp.parallel()) if (!receiveCpp.parallel())
@ -1083,7 +442,7 @@ void Foam::particle::hitCyclicAMIPatch
( (
receiveCpp.forwardT().size() == 1 receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0] ? receiveCpp.forwardT()[0]
: receiveCpp.forwardT()[patchFacei] : receiveCpp.forwardT()[receiveFacei]
); );
transformProperties(T); transformProperties(T);
} }
@ -1093,7 +452,7 @@ void Foam::particle::hitCyclicAMIPatch
( (
(receiveCpp.separation().size() == 1) (receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0] ? receiveCpp.separation()[0]
: receiveCpp.separation()[patchFacei] : receiveCpp.separation()[receiveFacei]
); );
transformProperties(-s); transformProperties(-s);
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -61,26 +61,26 @@ public:
passiveParticle passiveParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
particle(mesh, position, celli, tetFacei, tetPti) particle(mesh, coordinates, celli, tetFacei, tetPti)
{} {}
//- Construct from components, with searching for tetFace and
// tetPt unless disabled by doCellFacePt = false. //- Construct from a position and a cell, searching for the rest of the
// required topology
passiveParticle passiveParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli, const label celli
bool doCellFacePt = true
) )
: :
particle(mesh, position, celli, doCellFacePt) particle(mesh, position, celli)
{} {}
//- Construct from Istream //- Construct from Istream

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -321,7 +321,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
: -1 : -1
), ),
cellOccupancyPtr_(), cellOccupancyPtr_(),
cellLengthScale_(cbrt(mesh_.V())), cellLengthScale_(mag(cbrt(mesh_.V()))),
rho_(rho), rho_(rho),
U_(U), U_(U),
mu_(mu), mu_(mu),
@ -848,18 +848,14 @@ void Foam::KinematicCloud<CloudType>::updateMesh()
{ {
updateCellOccupancy(); updateCellOccupancy();
injectors_.updateMesh(); injectors_.updateMesh();
cellLengthScale_ = cbrt(mesh_.V()); cellLengthScale_ = mag(cbrt(mesh_.V()));
} }
template<class CloudType> template<class CloudType>
void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper) void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{ {
typedef typename particle::TrackingData<KinematicCloud<CloudType>> tdType; Cloud<parcelType>::autoMap(mapper);
tdType td(*this);
Cloud<parcelType>::template autoMap<tdType>(td, mapper);
updateMesh(); updateMesh();
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -333,11 +333,7 @@ void Foam::ReactingCloud<CloudType>::evolve()
template<class CloudType> template<class CloudType>
void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper) void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{ {
typedef typename particle::TrackingData<ReactingCloud<CloudType>> tdType; Cloud<parcelType>::autoMap(mapper);
tdType td(*this);
Cloud<parcelType>::template autoMap<tdType>(td, mapper);
this->updateMesh(); this->updateMesh();
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -259,12 +259,7 @@ void Foam::ReactingMultiphaseCloud<CloudType>::autoMap
const mapPolyMesh& mapper const mapPolyMesh& mapper
) )
{ {
typedef typename particle::TrackingData<ReactingMultiphaseCloud<CloudType>> Cloud<parcelType>::autoMap(mapper);
tdType;
tdType td(*this);
Cloud<parcelType>::template autoMap<tdType>(td, mapper);
this->updateMesh(); this->updateMesh();
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -483,11 +483,7 @@ void Foam::ThermoCloud<CloudType>::evolve()
template<class CloudType> template<class CloudType>
void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper) void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{ {
typedef typename particle::TrackingData<ThermoCloud<CloudType>> tdType; Cloud<parcelType>::autoMap(mapper);
tdType td(*this);
Cloud<parcelType>::template autoMap<tdType>(td, mapper);
this->updateMesh(); this->updateMesh();
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -163,22 +163,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline CollidingParcel inline CollidingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline CollidingParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline CollidingParcel inline CollidingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -63,13 +63,13 @@ template<class ParcelType>
inline Foam::CollidingParcel<ParcelType>::CollidingParcel inline Foam::CollidingParcel<ParcelType>::CollidingParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(owner, position, celli, tetFacei, tetPti), ParcelType(owner, coordinates, celli, tetFacei, tetPti),
f_(Zero), f_(Zero),
angularMomentum_(Zero), angularMomentum_(Zero),
torque_(Zero), torque_(Zero),
@ -82,6 +82,22 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(owner, position, celli),
f_(Zero),
angularMomentum_(Zero),
torque_(Zero),
collisionRecords_()
{}
template<class ParcelType>
inline Foam::CollidingParcel<ParcelType>::CollidingParcel
(
const polyMesh& owner,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -99,7 +115,7 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
ParcelType ParcelType
( (
owner, owner,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,

View File

@ -267,74 +267,50 @@ bool Foam::KinematicParcel<ParcelType>::move
const scalarField& cellLengthScale = td.cloud().cellLengthScale(); const scalarField& cellLengthScale = td.cloud().cellLengthScale();
const scalar maxCo = td.cloud().solution().maxCo(); const scalar maxCo = td.cloud().solution().maxCo();
scalar tEnd = (1.0 - p.stepFraction())*trackTime; while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
scalar dtMax = trackTime;
if (td.cloud().solution().transient())
{
dtMax *= maxCo;
}
bool tracking = true;
label nTrackingStalled = 0;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{ {
// Apply correction to position for reduced-D cases // Apply correction to position for reduced-D cases
meshTools::constrainToMeshCentre(mesh, p.position()); p.constrainToMeshCentre();
const point start(p.position()); // Cache the current position, cell and step-fraction
const point start = p.position();
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
// Cache the parcel current cell as this will change if a face is hit
const label celli = p.cell(); const label celli = p.cell();
const scalar sfrac = p.stepFraction();
const scalar magU = mag(U_); // Total displacement over the time-step
if (p.active() && tracking && (magU > ROOTVSMALL)) const vector s = trackTime*U_;
// Cell length scale
const scalar l = cellLengthScale[p.cell()];
// Fraction of the displacement to track in this loop. This is limited
// to ensure that the both the time and distance tracked is less than
// maxCo times the total value.
scalar f = 1 - p.stepFraction();
f = min(f, maxCo);
f = min(f, maxCo*l/max(SMALL*l, mag(s)));
if (p.active())
{ {
const scalar d = dt*magU; // Track to the next face
const scalar dCorr = min(d, maxCo*cellLengthScale[celli]); p.trackToFace(f*s, f, td);
dt *=
dCorr/d
*p.trackToFace(p.position() + dCorr*U_/magU, td);
} }
else
tEnd -= dt;
scalar newStepFraction = 1.0 - tEnd/trackTime;
if (tracking)
{ {
if // Abandon the track, and move to the end of the sub-step. If the
( // the mesh is moving, this will implicitly move the parcel.
mag(p.stepFraction() - newStepFraction) if (mesh.moving())
< particle::minStepFractionTol
)
{ {
nTrackingStalled++; WarningInFunction
<< "Tracking was abandoned on a moving mesh. Parcels may "
if (nTrackingStalled > maxTrackAttempts) << "move unphysically as a result." << endl;
{
tracking = false;
}
}
else
{
nTrackingStalled = 0;
} }
p.stepFraction() += f;
} }
p.stepFraction() = newStepFraction; const scalar dt = (p.stepFraction() - sfrac)*trackTime;
bool calcParcel = true;
if (!tracking && td.cloud().solution().steadyState())
{
calcParcel = false;
}
// Avoid problems with extremely small timesteps // Avoid problems with extremely small timesteps
if ((dt > ROOTVSMALL) && calcParcel) if (!td.cloud().solution().steadyState() && dt > ROOTVSMALL)
{ {
// Update cell based properties // Update cell based properties
p.setCellValues(td, dt, celli); p.setCellValues(td, dt, celli);
@ -347,7 +323,7 @@ bool Foam::KinematicParcel<ParcelType>::move
p.calc(td, dt, celli); p.calc(td, dt, celli);
} }
if (p.onBoundary() && td.keepParticle) if (p.onBoundaryFace() && td.keepParticle)
{ {
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())])) if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{ {

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -320,22 +320,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline KinematicParcel inline KinematicParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline KinematicParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline KinematicParcel inline KinematicParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -475,9 +484,6 @@ public:
// Helper functions // Helper functions
//- Return the index of the face used in the interpolation routine
inline label faceInterpolation() const;
//- Cell owner mass //- Cell owner mass
inline scalar massCell(const label celli) const; inline scalar massCell(const label celli) const;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -73,13 +73,13 @@ template<class ParcelType>
inline Foam::KinematicParcel<ParcelType>::KinematicParcel inline Foam::KinematicParcel<ParcelType>::KinematicParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(owner, position, celli, tetFacei, tetPti), ParcelType(owner, coordinates, celli, tetFacei, tetPti),
active_(true), active_(true),
typeId_(-1), typeId_(-1),
nParticle_(0), nParticle_(0),
@ -101,6 +101,31 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(owner, position, celli),
active_(true),
typeId_(-1),
nParticle_(0),
d_(0.0),
dTarget_(0.0),
U_(Zero),
rho_(0.0),
age_(0.0),
tTurb_(0.0),
UTurb_(Zero),
rhoc_(0.0),
Uc_(Zero),
muc_(0.0)
{}
template<class ParcelType>
inline Foam::KinematicParcel<ParcelType>::KinematicParcel
(
const polyMesh& owner,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -112,7 +137,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
const constantProperties& constProps const constantProperties& constProps
) )
: :
ParcelType(owner, position, celli, tetFacei, tetPti), ParcelType(owner, coordinates, celli, tetFacei, tetPti),
active_(true), active_(true),
typeId_(typeId), typeId_(typeId),
nParticle_(nParticle0), nParticle_(nParticle0),
@ -334,21 +359,6 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
} }
template<class ParcelType>
inline Foam::label Foam::KinematicParcel<ParcelType>::faceInterpolation() const
{
// Use volume-based interpolation if dealing with external faces
if (this->cloud().internalFace(this->face()))
{
return this->face();
}
else
{
return -1;
}
}
template<class ParcelType> template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell
( (

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -181,22 +181,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline MPPICParcel inline MPPICParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline MPPICParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline MPPICParcel inline MPPICParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -29,13 +29,13 @@ template<class ParcelType>
inline Foam::MPPICParcel<ParcelType>::MPPICParcel inline Foam::MPPICParcel<ParcelType>::MPPICParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(owner, position, celli, tetFacei, tetPti), ParcelType(owner, coordinates, celli, tetFacei, tetPti),
UCorrect_(Zero) UCorrect_(Zero)
{} {}
@ -45,6 +45,19 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
( (
const polyMesh& owner, const polyMesh& owner,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(owner, position, celli),
UCorrect_(Zero)
{}
template<class ParcelType>
inline Foam::MPPICParcel<ParcelType>::MPPICParcel
(
const polyMesh& owner,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -60,7 +73,7 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
ParcelType ParcelType
( (
owner, owner,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -270,23 +270,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, position and topology
// Other properties initialised as null // Other properties initialised as null
inline ReactingMultiphaseParcel inline ReactingMultiphaseParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline ReactingMultiphaseParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline ReactingMultiphaseParcel inline ReactingMultiphaseParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -68,13 +68,13 @@ template<class ParcelType>
inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(mesh, position, celli, tetFacei, tetPti), ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
YGas_(0), YGas_(0),
YLiquid_(0), YLiquid_(0),
YSolid_(0), YSolid_(0),
@ -87,6 +87,22 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(mesh, position, celli),
YGas_(0),
YLiquid_(0),
YSolid_(0),
canCombust_(0)
{}
template<class ParcelType>
inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
(
const polyMesh& mesh,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -108,7 +124,7 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
ParcelType ParcelType
( (
mesh, mesh,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -224,22 +224,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline ReactingParcel inline ReactingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline ReactingParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline ReactingParcel inline ReactingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -63,13 +63,13 @@ template<class ParcelType>
inline Foam::ReactingParcel<ParcelType>::ReactingParcel inline Foam::ReactingParcel<ParcelType>::ReactingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(mesh, position, celli, tetFacei, tetPti), ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
mass0_(0.0), mass0_(0.0),
Y_(0), Y_(0),
pc_(0.0) pc_(0.0)
@ -81,6 +81,21 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(mesh, position, celli),
mass0_(0.0),
Y_(0),
pc_(0.0)
{}
template<class ParcelType>
inline Foam::ReactingParcel<ParcelType>::ReactingParcel
(
const polyMesh& mesh,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -99,7 +114,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
ParcelType ParcelType
( (
mesh, mesh,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -280,22 +280,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline ThermoParcel inline ThermoParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline ThermoParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline ThermoParcel inline ThermoParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -74,13 +74,13 @@ template<class ParcelType>
inline Foam::ThermoParcel<ParcelType>::ThermoParcel inline Foam::ThermoParcel<ParcelType>::ThermoParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(mesh, position, celli, tetFacei, tetPti), ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
T_(0.0), T_(0.0),
Cp_(0.0), Cp_(0.0),
Tc_(0.0), Tc_(0.0),
@ -93,6 +93,22 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(mesh, position, celli),
T_(0.0),
Cp_(0.0),
Tc_(0.0),
Cpc_(0.0)
{}
template<class ParcelType>
inline Foam::ThermoParcel<ParcelType>::ThermoParcel
(
const polyMesh& mesh,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -110,7 +126,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
ParcelType ParcelType
( (
mesh, mesh,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -473,8 +473,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
meshTools::constrainToMeshCentre(mesh, pos); meshTools::constrainToMeshCentre(mesh, pos);
// Create a new parcel // Create a new parcel
parcelType* pPtr = parcelType* pPtr = new parcelType(mesh, pos, celli);
new parcelType(mesh, pos, celli, tetFacei, tetPti);
// Check/set new parcel thermo properties // Check/set new parcel thermo properties
cloud.setParcelThermoProperties(*pPtr, dt); cloud.setParcelThermoProperties(*pPtr, dt);
@ -590,8 +589,7 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
meshTools::constrainToMeshCentre(mesh, pos); meshTools::constrainToMeshCentre(mesh, pos);
// Create a new parcel // Create a new parcel
parcelType* pPtr = parcelType* pPtr = new parcelType(mesh, pos, celli);
new parcelType(mesh, pos, celli, tetFacei, tetPti);
// Check/set new parcel thermo properties // Check/set new parcel thermo properties
cloud.setParcelThermoProperties(*pPtr, 0.0); cloud.setParcelThermoProperties(*pPtr, 0.0);

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -146,16 +146,6 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
{ {
const label celli = injectorCellsPatch[j]; const label celli = injectorCellsPatch[j];
// The position could bein any tet of the decomposed cell,
// so arbitrarily choose the first face of the cell as the
// tetFace and the first point on the face after the base
// point as the tetPt. The tracking will pick the cell
// consistent with the motion in the first tracking step.
const label tetFacei = this->owner().mesh().cells()[celli][0];
const label tetPti = 1;
// const point& pos = this->owner().mesh().C()[celli];
const scalar offset = const scalar offset =
max max
( (
@ -166,14 +156,7 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
// Create a new parcel // Create a new parcel
parcelType* pPtr = parcelType* pPtr =
new parcelType new parcelType(this->owner().pMesh(), pos, celli);
(
this->owner().pMesh(),
pos,
celli,
tetFacei,
tetPti
);
// Check/set new parcel thermo properties // Check/set new parcel thermo properties
td.cloud().setParcelThermoProperties(*pPtr, 0.0); td.cloud().setParcelThermoProperties(*pPtr, 0.0);

View File

@ -449,7 +449,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
} }
// Perturb new parcels towards the owner cell centre // Perturb new parcels towards the owner cell centre
pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf); pPtr->track(0.5*rndGen_.sample01<scalar>()*(posC - posCf), 0);
pPtr->nParticle() = npNew[i]; pPtr->nParticle() = npNew[i];

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -85,18 +85,10 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
{ {
// Leapfrog tracking part // Leapfrog tracking part
scalar tEnd = (1.0 - stepFraction())*trackTime; while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{ {
// set the lagrangian time-step const scalar f = 1 - stepFraction();
scalar dt = min(dtMax, tEnd); trackToFace(f*trackTime*v_, f, td);
dt *= trackToFace(position() + dt*v_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
} }
} }
else if (td.part() == 2) else if (td.part() == 2)
@ -205,7 +197,7 @@ void Foam::molecule::transformProperties(const tensor& T)
rf_ = transform(T, rf_); rf_ = transform(T, rf_);
sitePositions_ = position_ + (T & (sitePositions_ - position_)); sitePositions_ = position() + (T & (sitePositions_ - position()));
siteForces_ = T & siteForces_; siteForces_ = T & siteForces_;
} }
@ -226,7 +218,7 @@ void Foam::molecule::transformProperties(const vector& separation)
void Foam::molecule::setSitePositions(const constantProperties& constProps) void Foam::molecule::setSitePositions(const constantProperties& constProps)
{ {
sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions()); sitePositions_ = position() + (Q_ & constProps.siteReferencePositions());
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -239,7 +239,7 @@ public:
inline molecule inline molecule
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -254,6 +254,24 @@ public:
const label id const label id
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology
inline molecule
(
const polyMesh& mesh,
const vector& position,
const label celli,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
);
//- Construct from Istream //- Construct from Istream
molecule molecule
( (

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -220,7 +220,7 @@ inline Foam::molecule::constantProperties::constantProperties
inline Foam::molecule::molecule inline Foam::molecule::molecule
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -236,7 +236,42 @@ inline Foam::molecule::molecule
) )
: :
particle(mesh, position, celli, tetFacei, tetPti), particle(mesh, coordinates, celli, tetFacei, tetPti),
Q_(Q),
v_(v),
a_(a),
pi_(pi),
tau_(tau),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(Zero),
special_(special),
id_(id),
siteForces_(constProps.nSites(), Zero),
sitePositions_(constProps.nSites())
{
setSitePositions(constProps);
}
inline Foam::molecule::molecule
(
const polyMesh& mesh,
const vector& position,
const label celli,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
)
:
particle(mesh, position, celli),
Q_(Q), Q_(Q),
v_(v), v_(v),
a_(a), a_(a),

View File

@ -746,17 +746,8 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition globalPosition
); );
label cell = -1; const label cell =
label tetFace = -1; mesh_.cellTree().findInside(globalPosition);
label tetPt = -1;
mesh_.findCellFacePt
(
globalPosition,
cell,
tetFace,
tetPt
);
if (findIndex(zone, cell) != -1) if (findIndex(zone, cell) != -1)
{ {
@ -764,8 +755,6 @@ void Foam::moleculeCloud::initialiseMolecules
( (
globalPosition, globalPosition,
cell, cell,
tetFace,
tetPt,
id, id,
tethered, tethered,
temperature, temperature,
@ -834,17 +823,11 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition globalPosition
); );
label cell = -1; const label cell =
label tetFace = -1; mesh_.cellTree().findInside
label tetPt = -1; (
globalPosition
mesh_.findCellFacePt );
(
globalPosition,
cell,
tetFace,
tetPt
);
if (findIndex(zone, cell) != -1) if (findIndex(zone, cell) != -1)
{ {
@ -852,8 +835,6 @@ void Foam::moleculeCloud::initialiseMolecules
( (
globalPosition, globalPosition,
cell, cell,
tetFace,
tetPt,
id, id,
tethered, tethered,
temperature, temperature,
@ -913,17 +894,11 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition globalPosition
); );
label cell = -1; const label cell =
label tetFace = -1; mesh_.cellTree().findInside
label tetPt = -1; (
globalPosition
mesh_.findCellFacePt );
(
globalPosition,
cell,
tetFace,
tetPt
);
if (findIndex(zone, cell) != -1) if (findIndex(zone, cell) != -1)
{ {
@ -931,8 +906,6 @@ void Foam::moleculeCloud::initialiseMolecules
( (
globalPosition, globalPosition,
cell, cell,
tetFace,
tetPt,
id, id,
tethered, tethered,
temperature, temperature,
@ -996,26 +969,12 @@ void Foam::moleculeCloud::createMolecule
( (
const point& position, const point& position,
label cell, label cell,
label tetFace,
label tetPt,
label id, label id,
bool tethered, bool tethered,
scalar temperature, scalar temperature,
const vector& bulkVelocity const vector& bulkVelocity
) )
{ {
if (cell == -1)
{
mesh_.findCellFacePt(position, cell, tetFace, tetPt);
}
if (cell == -1)
{
FatalErrorInFunction
<< "Position specified does not correspond to a mesh cell." << nl
<< abort(FatalError);
}
point specialPosition(Zero); point specialPosition(Zero);
label special = 0; label special = 0;
@ -1068,8 +1027,6 @@ void Foam::moleculeCloud::createMolecule
mesh_, mesh_,
position, position,
cell, cell,
tetFace,
tetPt,
Q, Q,
v, v,
Zero, Zero,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -114,8 +114,6 @@ private:
( (
const point& position, const point& position,
label cell, label cell,
label tetFace,
label tetPt,
label id, label id,
bool tethered, bool tethered,
scalar temperature, scalar temperature,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -43,34 +43,27 @@ bool Foam::solidParticle::move
td.switchProcessor = false; td.switchProcessor = false;
td.keepParticle = true; td.keepParticle = true;
const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh(); const polyBoundaryMesh& pbMesh = mesh().boundaryMesh();
scalar tEnd = (1.0 - stepFraction())*trackTime; while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
{ {
if (debug) if (debug)
{ {
Info<< "Time = " << mesh_.time().timeName() Info<< "Time = " << mesh().time().timeName()
<< " trackTime = " << trackTime << " trackTime = " << trackTime
<< " tEnd = " << tEnd
<< " steptFraction() = " << stepFraction() << endl; << " steptFraction() = " << stepFraction() << endl;
} }
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
// remember which cell the parcel is in const label celli = cell();
// since this will change if a face is hit const scalar sfrac = stepFraction();
label celli = cell();
dt *= trackToFace(position() + dt*U_, td); const scalar f = 1 - stepFraction();
trackToFace(f*trackTime*U_, f, td);
tEnd -= dt; const scalar dt = (stepFraction() - sfrac)*trackTime;
stepFraction() = 1.0 - tEnd/trackTime;
cellPointWeight cpw(mesh_, position(), celli, face()); cellPointWeight cpw(mesh(), position(), celli, face());
scalar rhoc = td.rhoInterp().interpolate(cpw); scalar rhoc = td.rhoInterp().interpolate(cpw);
vector Uc = td.UInterp().interpolate(cpw); vector Uc = td.UInterp().interpolate(cpw);
scalar nuc = td.nuInterp().interpolate(cpw); scalar nuc = td.nuInterp().interpolate(cpw);
@ -90,7 +83,7 @@ bool Foam::solidParticle::move
U_ = (U_ + dt*(Dc*Uc + (1.0 - rhoc/rhop)*td.g()))/(1.0 + dt*Dc); U_ = (U_ + dt*(Dc*Uc + (1.0 - rhoc/rhop)*td.g()))/(1.0 + dt*Dc);
if (onBoundary() && td.keepParticle) if (onBoundaryFace() && td.keepParticle)
{ {
if (isA<processorPolyPatch>(pbMesh[patch(face())])) if (isA<processorPolyPatch>(pbMesh[patch(face())]))
{ {
@ -133,7 +126,7 @@ void Foam::solidParticle::hitWallPatch
const tetIndices& tetIs const tetIndices& tetIs
) )
{ {
vector nw = tetIs.faceTri(mesh_).normal(); vector nw = tetIs.faceTri(mesh()).normal();
nw /= mag(nw); nw /= mag(nw);
scalar Un = U_ & nw; scalar Un = U_ & nw;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -130,7 +130,7 @@ public:
inline solidParticle inline solidParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -45,7 +45,7 @@ inline Foam::solidParticle::trackingData::trackingData
inline Foam::solidParticle::solidParticle inline Foam::solidParticle::solidParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -53,7 +53,7 @@ inline Foam::solidParticle::solidParticle
const vector& U const vector& U
) )
: :
particle(mesh, position, celli, tetFacei, tetPti), particle(mesh, coordinates, celli, tetFacei, tetPti),
d_(d), d_(d),
U_(U) U_(U)
{} {}

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -184,22 +184,31 @@ public:
// Constructors // Constructors
//- Construct from owner, position, and cloud owner //- Construct from mesh, coordinates and topology
// Other properties initialised as null // Other properties initialised as null
inline SprayParcel inline SprayParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology. Other properties are initialised as null.
inline SprayParcel
(
const polyMesh& mesh,
const vector& position,
const label celli
);
//- Construct from components //- Construct from components
inline SprayParcel inline SprayParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -107,15 +107,15 @@ template<class ParcelType>
inline Foam::SprayParcel<ParcelType>::SprayParcel inline Foam::SprayParcel<ParcelType>::SprayParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti const label tetPti
) )
: :
ParcelType(mesh, position, celli, tetFacei, tetPti), ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
d0_(this->d()), d0_(this->d()),
position0_(position), position0_(this->position()),
sigma_(0.0), sigma_(0.0),
mu_(0.0), mu_(0.0),
liquidCore_(0.0), liquidCore_(0.0),
@ -135,6 +135,31 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const vector& position,
const label celli
)
:
ParcelType(mesh, position, celli),
d0_(this->d()),
position0_(this->position()),
sigma_(0.0),
mu_(0.0),
liquidCore_(0.0),
KHindex_(0.0),
y_(0.0),
yDot_(0.0),
tc_(0.0),
ms_(0.0),
injector_(1.0),
tMom_(GREAT),
user_(0.0)
{}
template<class ParcelType>
inline Foam::SprayParcel<ParcelType>::SprayParcel
(
const polyMesh& mesh,
const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPti, const label tetPti,
@ -162,7 +187,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
ParcelType ParcelType
( (
mesh, mesh,
position, coordinates,
celli, celli,
tetFacei, tetFacei,
tetPti, tetPti,
@ -178,7 +203,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
constProps constProps
), ),
d0_(d0), d0_(d0),
position0_(position), position0_(this->position()),
sigma_(constProps.sigma0()), sigma_(constProps.sigma0()),
mu_(constProps.mu0()), mu_(constProps.mu0()),
liquidCore_(liquidCore), liquidCore_(liquidCore),

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -350,18 +350,7 @@ void Foam::meshRefinement::markFeatureCellLevel
{ {
const point& keepPoint = keepPoints[i]; const point& keepPoint = keepPoints[i];
label celli = -1; const label celli = mesh_.cellTree().findInside(keepPoint);
label tetFacei = -1;
label tetPtI = -1;
// Force construction of search tree even if processor holds no
// cells
(void)mesh_.cellTree();
if (mesh_.nCells())
{
mesh_.findCellFacePt(keepPoint, celli, tetFacei, tetPtI);
}
if (celli != -1) if (celli != -1)
{ {
@ -405,8 +394,6 @@ void Foam::meshRefinement::markFeatureCellLevel
mesh_, mesh_,
keepPoint, keepPoint,
celli, celli,
tetFacei,
tetPtI,
featureMesh.points()[pointi], // endpos featureMesh.points()[pointi], // endpos
featureLevel, // level featureLevel, // level
featI, // featureMesh featI, // featureMesh
@ -450,8 +437,6 @@ void Foam::meshRefinement::markFeatureCellLevel
mesh_, mesh_,
keepPoint, keepPoint,
celli, celli,
tetFacei,
tetPtI,
featureMesh.points()[pointi], // endpos featureMesh.points()[pointi], // endpos
featureLevel, // level featureLevel, // level
featI, // featureMesh featI, // featureMesh
@ -545,6 +530,7 @@ void Foam::meshRefinement::markFeatureCellLevel
label otherPointi = e.otherVertex(pointi); label otherPointi = e.otherVertex(pointi);
trackedParticle* tp(new trackedParticle(startTp)); trackedParticle* tp(new trackedParticle(startTp));
tp->start() = tp->position();
tp->end() = featureMesh.points()[otherPointi]; tp->end() = featureMesh.points()[otherPointi];
tp->j() = otherPointi; tp->j() = otherPointi;
tp->k() = edgeI; tp->k() = edgeI;
@ -605,6 +591,7 @@ void Foam::meshRefinement::markFeatureCellLevel
const edge& e = featureMesh.edges()[edgeI]; const edge& e = featureMesh.edges()[edgeI];
label otherPointi = e.otherVertex(pointi); label otherPointi = e.otherVertex(pointi);
tp.start() = tp.position();
tp.end() = featureMesh.points()[otherPointi]; tp.end() = featureMesh.points()[otherPointi];
tp.j() = otherPointi; tp.j() = otherPointi;
tp.k() = edgeI; tp.k() = edgeI;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -31,7 +31,7 @@ License
Foam::trackedParticle::trackedParticle Foam::trackedParticle::trackedParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPtI, const label tetPtI,
@ -42,7 +42,30 @@ Foam::trackedParticle::trackedParticle
const label k const label k
) )
: :
particle(mesh, position, celli, tetFacei, tetPtI), particle(mesh, coordinates, celli, tetFacei, tetPtI),
start_(position()),
end_(end),
level_(level),
i_(i),
j_(j),
k_(k)
{}
Foam::trackedParticle::trackedParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const point& end,
const label level,
const label i,
const label j,
const label k
)
:
particle(mesh, position, celli),
start_(this->position()),
end_(end), end_(end),
level_(level), level_(level),
i_(i), i_(i),
@ -64,7 +87,7 @@ Foam::trackedParticle::trackedParticle
{ {
if (is.format() == IOstream::ASCII) if (is.format() == IOstream::ASCII)
{ {
is >> end_; is >> start_ >> end_;
level_ = readLabel(is); level_ = readLabel(is);
i_ = readLabel(is); i_ = readLabel(is);
j_ = readLabel(is); j_ = readLabel(is);
@ -74,8 +97,8 @@ Foam::trackedParticle::trackedParticle
{ {
is.read is.read
( (
reinterpret_cast<char*>(&end_), reinterpret_cast<char*>(&start_),
sizeof(end_) + sizeof(level_) sizeof(start_) + sizeof(end_) + sizeof(level_)
+ sizeof(i_) + sizeof(j_) + sizeof(k_) + sizeof(i_) + sizeof(j_) + sizeof(k_)
); );
} }
@ -101,9 +124,8 @@ bool Foam::trackedParticle::move
td.switchProcessor = false; td.switchProcessor = false;
scalar tEnd = (1.0 - stepFraction())*trackTime; scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar dtMax = tEnd;
if (tEnd <= SMALL && onBoundary()) if (tEnd <= SMALL && onBoundaryFace())
{ {
// This is a hack to handle particles reaching their endpoint // This is a hack to handle particles reaching their endpoint
// on a processor boundary. If the endpoint is on a processor face // on a processor boundary. If the endpoint is on a processor face
@ -116,18 +138,13 @@ bool Foam::trackedParticle::move
{ {
td.keepParticle = true; td.keepParticle = true;
while (td.keepParticle && !td.switchProcessor && tEnd > SMALL) while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
{ {
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
// mark visited cell with max level. // mark visited cell with max level.
td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_); td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_);
dt *= trackToFace(end_, td); const scalar f = 1 - stepFraction();
trackToFace(f*(end_ - start_), f, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
} }
} }
@ -253,6 +270,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
if (os.format() == IOstream::ASCII) if (os.format() == IOstream::ASCII)
{ {
os << static_cast<const particle&>(p) os << static_cast<const particle&>(p)
<< token::SPACE << p.start_
<< token::SPACE << p.end_ << token::SPACE << p.end_
<< token::SPACE << p.level_ << token::SPACE << p.level_
<< token::SPACE << p.i_ << token::SPACE << p.i_
@ -264,8 +282,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
os << static_cast<const particle&>(p); os << static_cast<const particle&>(p);
os.write os.write
( (
reinterpret_cast<const char*>(&p.end_), reinterpret_cast<const char*>(&p.start_),
sizeof(p.end_) + sizeof(p.level_) sizeof(p.start_) + sizeof(p.end_) + sizeof(p.level_)
+ sizeof(p.i_) + sizeof(p.j_) + sizeof(p.k_) + sizeof(p.i_) + sizeof(p.j_) + sizeof(p.k_)
); );
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -64,6 +64,9 @@ class trackedParticle
{ {
// Private data // Private data
//- Start point to track from
point start_;
//- End point to track to //- End point to track to
point end_; point end_;
@ -120,7 +123,7 @@ public:
trackedParticle trackedParticle
( (
const polyMesh& mesh, const polyMesh& mesh,
const vector& position, const barycentric& coordinates,
const label celli, const label celli,
const label tetFacei, const label tetFacei,
const label tetPtI, const label tetPtI,
@ -131,6 +134,20 @@ public:
const label k const label k
); );
//- Construct from a position and a cell, searching for the rest of the
// required topology
trackedParticle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const point& end,
const label level,
const label i,
const label j,
const label k
);
//- Construct from Istream //- Construct from Istream
trackedParticle trackedParticle
( (
@ -170,6 +187,12 @@ public:
// Member Functions // Member Functions
//- Point to track from
point& start()
{
return start_;
}
//- Point to track to //- Point to track to
point& end() point& end()
{ {

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -48,9 +48,7 @@ void Foam::reconstructLagrangianPositions
forAll(meshes, i) forAll(meshes, i)
{ {
const labelList& cellMap = cellProcAddressing[i]; const labelList& cellMap = cellProcAddressing[i];
const labelList& faceMap = faceProcAddressing[i];
// faceProcAddressing not required currently.
// const labelList& faceMap = faceProcAddressing[i];
Cloud<passiveParticle> lpi(meshes[i], cloudName, false); Cloud<passiveParticle> lpi(meshes[i], cloudName, false);
@ -58,18 +56,21 @@ void Foam::reconstructLagrangianPositions
{ {
const passiveParticle& ppi = iter(); const passiveParticle& ppi = iter();
// // Inverting sign if necessary and subtracting 1 from const label mappedCell = cellMap[ppi.cell()];
// // faceProcAddressing
// label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1; // Inverting sign if necessary and subtracting 1 from
// faceProcAddressing
label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
lagrangianPositions.append lagrangianPositions.append
( (
new passiveParticle new passiveParticle
( (
mesh, mesh,
ppi.position(), ppi.coordinates(),
cellMap[ppi.cell()], mappedCell,
false mappedTetFace,
ppi.procTetPt(mesh, mappedCell, mappedTetFace)
) )
); );
} }

View File

@ -56,13 +56,15 @@ bool Foam::faceOnlySet::trackToBoundary
{ {
particle::TrackingData<passiveParticleCloud> trackData(particleCloud); particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
const point& trackPt = singleParticle.position(); point trackPt = singleParticle.position();
while(true) while(true)
{ {
point oldPoint = trackPt; point oldPoint = trackPt;
singleParticle.trackToFace(end_, trackData); singleParticle.trackToFace(end_ - start_, 0, trackData);
trackPt = singleParticle.position();
if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist) if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
{ {
@ -78,7 +80,7 @@ bool Foam::faceOnlySet::trackToBoundary
// End reached // End reached
return false; return false;
} }
else if (singleParticle.onBoundary()) else if (singleParticle.onBoundaryFace())
{ {
// Boundary reached // Boundary reached
return true; return true;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -56,29 +56,16 @@ bool Foam::polyLineSet::trackToBoundary
{ {
particle::TrackingData<passiveParticleCloud> trackData(particleCloud); particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
// Alias
const point& trackPt = singleParticle.position();
while (true) while (true)
{ {
// Local geometry info // Local geometry info
const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI]; const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
const scalar smallDist = mag(tol*offset); const scalar smallDist = mag(tol*offset);
point oldPos = trackPt; singleParticle.track(offset, 0);
label facei = -1; const point trackPt = singleParticle.position();
do
{
singleParticle.stepFraction() = 0;
singleParticle.track(sampleCoords_[sampleI+1], trackData);
}
while
(
!singleParticle.onBoundary()
&& (mag(trackPt - oldPos) < smallDist)
);
if (singleParticle.onBoundary()) if (singleParticle.onBoundaryFace())
{ {
//Info<< "trackToBoundary : reached boundary" //Info<< "trackToBoundary : reached boundary"
// << " trackPt:" << trackPt << endl; // << " trackPt:" << trackPt << endl;
@ -94,7 +81,7 @@ bool Foam::polyLineSet::trackToBoundary
// << endl; // << endl;
samplingPts.append(trackPt); samplingPts.append(trackPt);
samplingCells.append(singleParticle.cell()); samplingCells.append(singleParticle.cell());
samplingFaces.append(facei); samplingFaces.append(singleParticle.face());
// trackPt is at sampleI+1 // trackPt is at sampleI+1
samplingCurveDist.append(1.0*(sampleI+1)); samplingCurveDist.append(1.0*(sampleI+1));

View File

@ -209,17 +209,25 @@ Foam::point Foam::sampledSet::pushIn
label tetPtI; label tetPtI;
mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI); mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
// This is the tolerance that was defined as a static constant of the
// particle class. It is no longer used by particle, following the switch to
// barycentric tracking. The only place in which the tolerance is now used
// is here. I'm not sure what the purpose of this code is, but it probably
// wants removing. It is doing tet-searches for a particle position. This
// should almost certainly be left to the particle class.
const scalar trackingCorrectionTol = 1e-5;
if (tetFacei == -1 || tetPtI == -1) if (tetFacei == -1 || tetPtI == -1)
{ {
newPosition = facePt; newPosition = facePt;
label trap(1.0/particle::trackingCorrectionTol + 1); label trap(1.0/trackingCorrectionTol + 1);
label iterNo = 0; label iterNo = 0;
do do
{ {
newPosition += particle::trackingCorrectionTol*(cC - facePt); newPosition += trackingCorrectionTol*(cC - facePt);
mesh().findTetFacePt mesh().findTetFacePt
( (

View File

@ -94,8 +94,7 @@ bool Foam::uniformSet::trackToBoundary
const vector smallVec = tol*offset; const vector smallVec = tol*offset;
const scalar smallDist = mag(smallVec); const scalar smallDist = mag(smallVec);
// Alias point trackPt = singleParticle.position();
const point& trackPt = singleParticle.position();
particle::TrackingData<passiveParticleCloud> trackData(particleCloud); particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
@ -153,32 +152,10 @@ bool Foam::uniformSet::trackToBoundary
<< " to:" << samplePt << endl; << " to:" << samplePt << endl;
} }
point oldPos = trackPt; singleParticle.track(samplePt - trackPt, 0);
label facei = -1; trackPt = singleParticle.position();
do
{
singleParticle.stepFraction() = 0;
singleParticle.track(samplePt, trackData);
if (debug) if (singleParticle.onBoundaryFace())
{
Pout<< "Result of tracking "
<< " trackPt:" << trackPt
<< " trackCelli:" << singleParticle.cell()
<< " trackFacei:" << singleParticle.face()
<< " onBoundary:" << singleParticle.onBoundary()
<< " samplePt:" << samplePt
<< " smallDist:" << smallDist
<< endl;
}
}
while
(
!singleParticle.onBoundary()
&& (mag(trackPt - oldPos) < smallDist)
);
if (singleParticle.onBoundary())
{ {
//Pout<< "trackToBoundary : reached boundary" << endl; //Pout<< "trackToBoundary : reached boundary" << endl;
if (mag(trackPt - samplePt) < smallDist) if (mag(trackPt - samplePt) < smallDist)
@ -188,7 +165,7 @@ bool Foam::uniformSet::trackToBoundary
// Reached samplePt on boundary // Reached samplePt on boundary
samplingPts.append(trackPt); samplingPts.append(trackPt);
samplingCells.append(singleParticle.cell()); samplingCells.append(singleParticle.cell());
samplingFaces.append(facei); samplingFaces.append(singleParticle.face());
samplingCurveDist.append(mag(trackPt - start_)); samplingCurveDist.append(mag(trackPt - start_));
} }

View File

@ -47,7 +47,6 @@ runTimeModifiable true;
functions functions
{ {
#include "streamLines" #include "streamLines"
#include "wallBoundedStreamLines"
#include "cuttingPlane" #include "cuttingPlane"
#include "forceCoeffs" #include "forceCoeffs"
} }

View File

@ -1,100 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
// Interpolate U to create near-wall UNear
near
{
// Where to load it from
libs ("libfieldFunctionObjects.so");
type nearWallFields;
// Output every
writeControl writeTime;
// Calculate every
executeControl writeTime;
// Fields to be sampled. Per field original name and mapped field to
// create.
// Note: fields only get updated when writing!
fields
(
(U UNear)
);
// Patches/groups to sample (regular expressions)
patches (motorBikeGroup);
// Distance to sample
distance 0.001;
}
// Use UNear to track along wall
wallBoundedStreamLines
{
// Where to load it from (if not already in solver)
libs ("libfieldFunctionObjects.so");
type wallBoundedStreamLine;
// Output every
writeControl writeTime;
setFormat vtk; //gnuplot; //xmgr; //raw; //jplot;
// Velocity field to use for tracking.
U UNear;
// Interpolation method. Default is cellPoint.
// interpolationScheme pointMVC;
// Tracked forwards (+U) or backwards (-U)
trackForward true;
interpolationScheme cellPoint;
// Names of fields to sample. Should contain above velocity field!
fields (p U k UNear);
// Steps particles can travel before being removed
lifeTime 100;
// Cloud name to use
cloudName wallBoundedParticleTracks;
// Seeding method.
seedSampleSet
{
type patchSeed;
uniformCoeffs
{
axis x; //distance;
start (0.0035 0.0999 0.0001);
end (0.0035 0.0999 0.0099);
nPoints 20;
}
cloudCoeffs
{
axis x; //distance;
points ((0.351516548679288 -0.0116085375585099 1.24));
}
patchSeedCoeffs
{
patches (motorBikeGroup);
axis x; //distance;
maxPoints 20000;
}
}
}
// ************************************************************************* //