diff --git a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
index aec4c643045d6b8303571eeee84531d17622b1d0..ded2af139beaed6df96de9a1e8de6d76b86db7b9 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
@@ -77,6 +77,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
+ // Store the particle positions
+ kinematicCloud.storeGlobalPositions();
+
mesh.update();
// Calculate absolute flux from the mapped surface velocity
diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
index be0905796ba3c19f8b98aefd7ddcc8802cfb7cc1..00c66823eafde1ddd85295da714ac2acd0b6721c 100644
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -68,6 +68,8 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
+ kinematicCloud.storeGlobalPositions();
+
mesh.update();
U.correctBoundaryConditions();
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
index 16d2d524c5b5ebbd5a2e8fc7b1be7c591e39affe..40acccae9e10d6fcda7a9891163543792a335f4b 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -91,6 +91,9 @@ int main(int argc, char *argv[])
// Store momentum to set rhoUf for introduced faces.
volVectorField rhoU("rhoU", rho*U);
+ // Store the particle positions
+ parcels.storeGlobalPositions();
+
// Do any mesh changes
mesh.update();
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..430b1df24fc54babe6e16dcaa62e81cba8d37da3
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files
@@ -0,0 +1,3 @@
+uncoupledKinematicParcelDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelDyMFoam
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..2d385f9853cb5a2bd407b6ac79fc743395f9334f
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@@ -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
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..18f8e638693e0aafe48347a90b8f8f694e013ebd
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
index efb7dd2a235dff43deb230588bd961c867862632..98457dc2b4e5e8ac5b80a94dd875f4aa4ff9a961 100644
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
+++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -47,13 +47,12 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
{
label pi = 0;
- // faceProcAddressing not required currently
- // labelList decodedProcFaceAddressing(faceProcAddressing.size());
+ labelList decodedProcFaceAddressing(faceProcAddressing.size());
- // forAll(faceProcAddressing, i)
- // {
- // decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
- // }
+ forAll(faceProcAddressing, i)
+ {
+ decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
+ }
forAll(cellProcAddressing, procCelli)
{
@@ -68,27 +67,28 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
const indexedParticle& ppi = *iter();
particleIndices_[pi++] = ppi.index();
- // label mappedTetFace = findIndex
- // (
- // decodedProcFaceAddressing,
- // ppi.tetFace()
- // );
+ label mappedTetFace = findIndex
+ (
+ decodedProcFaceAddressing,
+ ppi.tetFace()
+ );
- // if (mappedTetFace == -1)
- // {
- // FatalErrorInFunction
- // << "Face lookup failure." << nl
- // << abort(FatalError);
- // }
+ if (mappedTetFace == -1)
+ {
+ FatalErrorInFunction
+ << "Face lookup failure." << nl
+ << abort(FatalError);
+ }
positions_.append
(
new passiveParticle
(
procMesh,
- ppi.position(),
+ ppi.coordinates(),
procCelli,
- false
+ mappedTetFace,
+ ppi.procTetPt(procMesh, procCelli, mappedTetFace)
)
);
}
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
index ad80ba3d07c632dc2d8ffdcdbb64159c3fa3724b..45cd8e4bb5bdfba633c6a5b1d34032536c5e631b 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -98,7 +98,6 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
- const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs
@@ -182,15 +181,17 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
new passiveParticle
(
meshTarget,
- targetCc[targetCells[i]],
- targetCells[i]
+ barycentric(1, 0, 0, 0),
+ targetCells[i],
+ meshTarget.cells()[targetCells[i]][0],
+ 1
)
);
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.
foundCell = true;
@@ -233,11 +234,16 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
{
unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI);
- iter().cell() = targetCell;
targetParcels.addParticle
(
- sourceParcels.remove(&iter())
+ new passiveParticle
+ (
+ meshTarget,
+ iter().position(),
+ targetCell
+ )
);
+ sourceParcels.remove(&iter());
}
}
sourceParticleI++;
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
index 06f7edc9106e562c7a8f4d0509396a7d243be52f..2f41a6eee60cebfed591f02c3f9f6016f2b2a7a1 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -89,8 +89,6 @@ void mapLagrangian(const meshToMesh& interp)
const polyMesh& meshTarget = interp.tgtRegion();
const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
- const pointField& targetCc = meshTarget.cellCentres();
-
fileNameList cloudDirs
(
readDir
@@ -173,15 +171,17 @@ void mapLagrangian(const meshToMesh& interp)
new passiveParticle
(
meshTarget,
- targetCc[targetCells[i]],
- targetCells[i]
+ barycentric(1, 0, 0, 0),
+ targetCells[i],
+ meshTarget.cells()[targetCells[i]][0],
+ 1
)
);
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.
foundCell = true;
@@ -224,11 +224,16 @@ void mapLagrangian(const meshToMesh& interp)
{
unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI);
- iter().cell() = targetCell;
targetParcels.addParticle
(
- sourceParcels.remove(&iter())
+ new passiveParticle
+ (
+ meshTarget,
+ iter().position(),
+ targetCell
+ )
);
+ sourceParcels.remove(&iter());
}
}
sourceParticleI++;
diff --git a/etc/caseDicts/meshQualityDict b/etc/caseDicts/meshQualityDict
index e69207563a931bfcd4962fb610e9cd63dcb661aa..fe5a1f321add8da1933c2ac1b50031c83ed0dcd9 100644
--- a/etc/caseDicts/meshQualityDict
+++ b/etc/caseDicts/meshQualityDict
@@ -33,8 +33,7 @@ minVol 1e-13;
//- Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
-// the cell centre. This has to be a positive number for tracking
-// to work. Set to very negative number (e.g. -1E30) to
+// the cell centre. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index d9f2104c1831186f3e8a24e2e052509245f60f68..f79161fa78c917cb415f88efcb08cef72fd0c700 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -133,6 +133,8 @@ $(spatialVectorAlgebra)/CompactSpatialTensor/compactSpatialTensor/compactSpatial
primitives/polynomialEqns/cubicEqn/cubicEqn.C
primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
+primitives/Barycentric/barycentric/barycentric.C
+
containers/HashTables/HashTable/HashTableCore.C
containers/HashTables/StaticHashTable/StaticHashTableCore.C
containers/Lists/SortableList/ParSortableListName.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 8437bcde5c0bab9431149e21d6773d8b9c140355..52b8c4c54c0c9d7f06b6e49d1eca697dc9c5c6f4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 * * * * * * * * * * * * * * //
Foam::polyMesh::polyMesh(const IOobject& io)
@@ -207,7 +230,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
comm_(UPstream::worldComm),
geometricD_(Zero),
solutionD_(Zero),
- tetBasePtIsPtr_(nullptr),
+ tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr),
pointZones_
(
@@ -401,7 +424,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm),
geometricD_(Zero),
solutionD_(Zero),
- tetBasePtIsPtr_(nullptr),
+ tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr),
pointZones_
(
@@ -552,7 +575,7 @@ Foam::polyMesh::polyMesh
comm_(UPstream::worldComm),
geometricD_(Zero),
solutionD_(Zero),
- tetBasePtIsPtr_(nullptr),
+ tetBasePtIsPtr_(readTetBasePtIs()),
cellTreePtr_(nullptr),
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())
{
@@ -828,8 +851,17 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
tetBasePtIsPtr_.reset
(
- new labelList
+ new labelIOList
(
+ IOobject
+ (
+ "tetBasePtIs",
+ instance(),
+ meshSubDir,
+ *this,
+ IOobject::READ_IF_PRESENT,
+ IOobject::NO_WRITE
+ ),
polyMeshTetDecomposition::findFaceBasePts(*this)
)
);
@@ -1105,6 +1137,13 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
points_.instance() = time().timeName();
points_.eventNo() = getEvent();
+ if (tetBasePtIsPtr_.valid())
+ {
+ tetBasePtIsPtr_().writeOpt() = IOobject::AUTO_WRITE;
+ tetBasePtIsPtr_().instance() = time().timeName();
+ tetBasePtIsPtr_().eventNo() = getEvent();
+ }
+
tmp<scalarField> sweptVols = primitiveMesh::movePoints
(
points_,
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 1f61f7382835afffbad0360ec77c2c3447e9be83..00ff5aa2f9ee85e8ef3522cd22acb1e24f6d7d2e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -150,7 +150,7 @@ private:
mutable Vector<label> solutionD_;
//- Base point for face decomposition into tets
- mutable autoPtr<labelList> tetBasePtIsPtr_;
+ mutable autoPtr<labelIOList> tetBasePtIsPtr_;
//- Search tree to allow spatial cell searching
mutable autoPtr<indexedOctree<treeDataCell>> cellTreePtr_;
@@ -208,6 +208,9 @@ private:
// polyhedral information
void calcCellShapes() const;
+ //- Read and return the tetBasePtIs
+ autoPtr<labelIOList> readTetBasePtIs() const;
+
// Helper functions for constructor from cell shapes
@@ -449,7 +452,7 @@ public:
label nSolutionD() const;
//- Return the tetBasePtIs
- const labelList& tetBasePtIs() const;
+ const labelIOList& tetBasePtIs() const;
//- Return the cell search tree
const indexedOctree<treeDataCell>& cellTree() const;
@@ -605,8 +608,8 @@ public:
//- Clear primitive data (points, faces and cells)
void clearPrimitives();
- //- Clear geometry not used for CFD (cellTree, tetBasePtIs)
- void clearAdditionalGeom();
+ //- Clear tet base points
+ void clearTetBasePtIs();
//- Clear cell tree data
void clearCellTree();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index fa1bdc8bdd10b566bb792533cb5b909436782da1..73c4537ef24ac8d96792b0b4ee9e332b2ba6b1a3 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -69,22 +69,6 @@ void Foam::polyMesh::clearGeom()
geometricD_ = 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
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()
{
if (debug)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
index 7803f094bee53edbf47d286bdcdc3d1916c7cd41..544ef0ecfc770b3f7a5fb8a9e5aa403c5f5af2c4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -59,6 +59,12 @@ void Foam::polyMesh::setInstance(const fileName& inst)
cellZones_.writeOpt() = IOobject::AUTO_WRITE;
cellZones_.instance() = inst;
+
+ if (tetBasePtIsPtr_.valid())
+ {
+ tetBasePtIsPtr_->writeOpt() = IOobject::AUTO_WRITE;
+ tetBasePtIsPtr_->instance() = inst;
+ }
}
@@ -387,6 +393,9 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
cellZones_.set(czI, newCellZones[czI].clone(cellZones_));
}
+ // Re-read tet base points
+ tetBasePtIsPtr_ = readTetBasePtIs();
+
if (boundaryChanged)
{
@@ -439,6 +448,13 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
points_.transfer(newPoints);
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.)
boundary_.calcGeometry();
diff --git a/src/OpenFOAM/primitives/Barycentric/Barycentric.H b/src/OpenFOAM/primitives/Barycentric/Barycentric.H
new file mode 100644
index 0000000000000000000000000000000000000000..e430401cf8e64d445e22d4ede211ebdfa8f2b8ef
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/Barycentric.H
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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/>.
+
+Class
+ Foam::Barycentric
+
+Description
+ Templated 3D Barycentric derived from VectorSpace. Has 4 components, one of
+ which is redundant.
+
+SourceFiles
+ BarycentricI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Barycentric_H
+#define Barycentric_H
+
+#include "VectorSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class Barycentric Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Cmpt>
+class Barycentric
+:
+ public VectorSpace<Barycentric<Cmpt>, Cmpt, 4>
+{
+public:
+
+ //- Equivalent type of labels used for valid component indexing
+ 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
+
+ //- Construct null
+ inline Barycentric();
+
+ //- Construct initialized to zero
+ inline Barycentric(const Foam::zero);
+
+ //- Construct given four components
+ inline Barycentric
+ (
+ const Cmpt& va,
+ const Cmpt& vb,
+ const Cmpt& vc,
+ const Cmpt& vd
+ );
+
+
+ // Member Functions
+
+ // Access
+
+ 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();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BarycentricI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricI.H b/src/OpenFOAM/primitives/Barycentric/BarycentricI.H
new file mode 100644
index 0000000000000000000000000000000000000000..ad649cdac5fbb34c9c0606d24b260ad08a9bd8e5
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricI.H
@@ -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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H b/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H
new file mode 100644
index 0000000000000000000000000000000000000000..2d83fafd8fbcfb56d3669d5aa1620cf31979d134
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H
@@ -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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H b/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H
new file mode 100644
index 0000000000000000000000000000000000000000..84e461377f72c1921f288e1705d77677eaa73fa5
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H
@@ -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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C
new file mode 100644
index 0000000000000000000000000000000000000000..f42aa1a9fe3f8384234c050e6bca3816862fca3c
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C
@@ -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>()
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H
new file mode 100644
index 0000000000000000000000000000000000000000..34cef96a1fedbdd0302331352d4b476d140a3ac5
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H
@@ -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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H b/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H
new file mode 100644
index 0000000000000000000000000000000000000000..3bd3f06195dee8012db7cd8139f1c04615c82d6e
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::barycentricTensor
+
+Description
+ A scalar version of the templated BarycentricTensor
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef barycentricTensor_H
+#define barycentricTensor_H
+
+#include "scalar.H"
+#include "BarycentricTensor.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef BarycentricTensor<scalar> barycentricTensor;
+
+
+template<>
+inline bool contiguous<barycentricTensor>()
+{
+ return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
index db923761452ff0235937a111af8b366b171fbc07..633fe479f7eff33da5990697d7e40a8b7db0057a 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
@@ -92,6 +92,9 @@ public:
//- Evaluate at x
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
inline scalar error(const scalar x) const;
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
index eb6a3283e45156413c9251bb6ac232f95e4766ff..da5fc3cba14ef44e5f226d1bcdaef220455d540a 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
@@ -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
{
return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c()));
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
index debc52cc01ddb7875d6ca9e07bf93154d8a05d84..7f88336e631b290349d6ff907e3bea90a2c5397f 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
@@ -81,6 +81,9 @@ public:
//- Evaluate at x
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
inline scalar error(const scalar x) const;
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
index e86ccbdd0a6d9105504a3552eea7e40c97a9f493..17bfdf452da1f4e0a576b199770a0b5fbe20ea86 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
@@ -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
{
return mag(SMALL*x*a());
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
index 1327dc15c943b8acb0429a2658d69bbe445a8a22..c676227de18acbc1eb361bf76ea2c199e1468993 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
@@ -84,6 +84,9 @@ public:
//- Evaluate at x
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
inline scalar error(const scalar x) const;
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
index 4b4aacae614a57de75db138ed46ead5b1fb23f4e..7672a74d72458e4b70d5b2912a79c249bb0a8035 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
@@ -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
{
return mag(SMALL*x*(x*2*a() + b()));
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index ee586e29690afa1c49d73def4f8a09f81bf7ce70..e7d6f558d19a09ba45c6c03f2cc52b396c38f174 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -727,7 +727,7 @@ void Foam::motionSmootherAlgo::movePoints()
{
// Make sure to clear out tetPtIs since used in checks (sometimes, should
// really check)
- mesh_.clearAdditionalGeom();
+ mesh_.clearTetBasePtIs();
pp_.movePoints(mesh_.points());
}
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index b44dc3819f9a943fc135f45021d6dfcf89b200c9..5161f46cf23f9c5ec1e091f41d6c2d24430519a7 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -25,11 +25,6 @@ streamLine/streamLineBase.C
streamLine/streamLineParticle.C
streamLine/streamLineParticleCloud.C
-wallBoundedStreamLine/wallBoundedStreamLine.C
-wallBoundedStreamLine/wallBoundedStreamLineParticle.C
-wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
-wallBoundedStreamLine/wallBoundedParticle.C
-
surfaceInterpolate/surfaceInterpolate.C
regionSizeDistribution/regionSizeDistribution.C
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C
index 6c69be4f18052c80d929e1d62f6584f7a1dfb8c7..e8efde11669a75ded2526bb394eae2c97e0c4826 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.C
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -30,7 +30,7 @@ License
Foam::findCellParticle::findCellParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
@@ -38,7 +38,24 @@ Foam::findCellParticle::findCellParticle
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),
data_(data)
{}
@@ -57,15 +74,15 @@ Foam::findCellParticle::findCellParticle
{
if (is.format() == IOstream::ASCII)
{
- is >> end_;
+ is >> start_ >> end_;
data_ = readLabel(is);
}
else
{
is.read
(
- reinterpret_cast<char*>(&end_),
- sizeof(end_) + sizeof(data_)
+ reinterpret_cast<char*>(&start_),
+ sizeof(start_) + sizeof(end_) + sizeof(data_)
);
}
}
@@ -85,21 +102,13 @@ bool Foam::findCellParticle::move
td.switchProcessor = false;
td.keepParticle = true;
- scalar tEnd = (1.0 - stepFraction())*maxTrackLen;
- scalar dtMax = tEnd;
-
- while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+ while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
{
- // set the lagrangian time-step
- scalar dt = min(dtMax, tEnd);
-
- dt *= trackToFace(end_, td);
-
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/maxTrackLen;
+ const scalar f = 1 - stepFraction();
+ trackToFace(f*(end_ - start_), f, td);
}
- if (tEnd < SMALL || !td.keepParticle)
+ if (stepFraction() == 1 || !td.keepParticle)
{
// Hit endpoint or patch. If patch hit could do fancy stuff but just
// to use the patch point is good enough for now.
@@ -209,6 +218,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
if (os.format() == IOstream::ASCII)
{
os << static_cast<const particle&>(p)
+ << token::SPACE << p.start_
<< token::SPACE << p.end_
<< token::SPACE << p.data_;
}
@@ -217,8 +227,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
os << static_cast<const particle&>(p);
os.write
(
- reinterpret_cast<const char*>(&p.end_),
- sizeof(p.end_) + sizeof(p.data_)
+ reinterpret_cast<const char*>(&p.start_),
+ sizeof(p.start_) + sizeof(p.end_) + sizeof(p.data_)
);
}
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.H b/src/functionObjects/field/nearWallFields/findCellParticle.H
index d4bb541807d153f17bf8974118a5588bcc57b0a2..ccc44e923880cae4aea2034c95529864f8b81d1c 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.H
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -63,6 +63,9 @@ class findCellParticle
{
// Private data
+ //- Start point to track from
+ point start_;
+
//- End point to track to
point end_;
@@ -119,7 +122,7 @@ public:
findCellParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
@@ -127,6 +130,17 @@ public:
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
findCellParticle
(
@@ -166,18 +180,42 @@ public:
// Member Functions
+ //- Point to track from
+ const point& start() const
+ {
+ return start_;
+ }
+
+ //- Point to track from
+ point& start()
+ {
+ return start_;
+ }
+
//- Point to track to
const point& end() const
{
return end_;
}
+ //- Point to track to
+ point& end()
+ {
+ return end_;
+ }
+
//- Transported label
label data() const
{
return data_;
}
+ //- Transported label
+ label& data()
+ {
+ return data_;
+ }
+
// Tracking
diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C
index abf2e0d4002921fde9a6bf7fa7132a6f08356f75..515f22763ab851221a80b1f1e4e3fcbb3c9229cb 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFields.C
+++ b/src/functionObjects/field/nearWallFields/nearWallFields.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -126,32 +126,19 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
const point end = start-distance_*nf[patchFacei];
- if (tetFacei == -1)
- {
- WarningInFunction << "Did not find point " << start
- << " inside cell " << celli
- << ". Not seeding particle originating from face centre "
- << mesh_.faceCentres()[meshFacei]
- << " with cell centre " << mesh_.cellCentres()[celli]
- << endl;
- }
- else
- {
- // Add to cloud. Add originating face as passive data
- cloud.addParticle
+ // Add a particle to the cloud with originating face as passive data
+ cloud.addParticle
+ (
+ new findCellParticle
(
- new findCellParticle
- (
- mesh_,
- start,
- celli,
- tetFacei,
- tetPti,
- end,
- globalWalls.toGlobal(nPatchFaces) // passive data
- )
- );
- }
+ mesh_,
+ start,
+ -1,
+ end,
+ globalWalls.toGlobal(nPatchFaces) // passive data
+ )
+ );
+
nPatchFaces++;
}
}
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.C b/src/functionObjects/field/streamLine/streamLineParticle.C
index 04189af71edf0f71e201a7d937198e04dc8372d1..c587a7279cafa0e6b579550bd2de8f302c9dce5a 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/functionObjects/field/streamLine/streamLineParticle.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -26,35 +26,8 @@ License
#include "streamLineParticle.H"
#include "vectorFieldIOField.H"
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-// defineParticleTypeNameAndDebug(streamLineParticle, 0);
-}
-
-
// * * * * * * * * * * * * * 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
(
const trackingData& td,
@@ -129,7 +102,6 @@ Foam::streamLineParticle::streamLineParticle
{
if (readFields)
{
- //if (is.format() == IOstream::ASCII)
List<scalarList> sampledScalars;
List<vectorList> sampledVectors;
@@ -169,31 +141,22 @@ Foam::streamLineParticle::streamLineParticle
bool Foam::streamLineParticle::move
(
trackingData& td,
- const scalar trackTime
+ const scalar
)
{
- streamLineParticle& p = static_cast<streamLineParticle&>(*this);
-
td.switchProcessor = false;
td.keepParticle = true;
- scalar tEnd = (1.0 - stepFraction())*trackTime;
- scalar maxDt = mesh_.bounds().mag();
+ const scalar maxDt = mesh().bounds().mag();
- while
- (
- td.keepParticle
- && !td.switchProcessor
- && lifeTime_ > 0
- )
+ while (td.keepParticle && !td.switchProcessor && lifeTime_ > 0)
{
- // set the lagrangian time-step
scalar dt = maxDt;
// Cross cell in steps:
// - at subiter 0 calculate dt to cross cell in nSubCycle steps
// - 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_;
@@ -219,37 +182,27 @@ bool Foam::streamLineParticle::move
if (td.trackLength_ < GREAT)
{
+ // No sub-cycling. Track a set length on each step.
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
- dt = calcSubCycleDeltaT(td, dt, U);
+ // Sub-cycling. Cross the cell in nSubCycle steps.
+ particle copy(*this);
+ copy.trackToFace(maxDt*U, 1);
+ dt *= (copy.stepFraction() - stepFraction())/td.nSubCycle_;
}
else if (subIter == td.nSubCycle_ - 1)
{
- // Do full step on last subcycle
+ // Sub-cycling. Track the whole cell on the last step.
dt = maxDt;
}
-
- scalar fraction = trackToFace(position() + dt*U, td);
- dt *= fraction;
-
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/trackTime;
-
- if (tEnd <= ROOTVSMALL)
- {
- // Force removal
- lifeTime_ = 0;
- }
+ trackToFace(dt*U, 0, td);
if
(
- face() != -1
+ onFace()
|| !td.keepParticle
|| td.switchProcessor
|| lifeTime_ == 0
@@ -260,17 +213,16 @@ bool Foam::streamLineParticle::move
}
}
-
if (!td.keepParticle || lifeTime_ == 0)
{
if (lifeTime_ == 0)
{
+ // Failure exit. Particle stagnated or it's life ran out.
if (debug)
{
- Pout<< "streamLineParticle : Removing stagnant particle:"
- << p.position()
- << " sampled positions:" << sampledPositions_.size()
- << endl;
+ Pout<< "streamLineParticle: Removing stagnant particle:"
+ << position() << " sampled positions:"
+ << sampledPositions_.size() << endl;
}
td.keepParticle = false;
}
@@ -282,29 +234,25 @@ bool Foam::streamLineParticle::move
if (debug)
{
- Pout<< "streamLineParticle : Removing particle:"
- << p.position()
+ Pout<< "streamLineParticle: Removing particle:" << 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]);
@@ -428,18 +376,11 @@ void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
);
c.checkFieldIOobject(c, sampledPositions);
-// vectorFieldIOField sampleVelocity
-// (
-// c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
-// );
-// c.checkFieldIOobject(c, sampleVelocity);
-
label i = 0;
forAllIter(Cloud<streamLineParticle>, c, iter)
{
iter().lifeTime_ = lifeTime[i];
iter().sampledPositions_.transfer(sampledPositions[i]);
-// iter().sampleVelocity_.transfer(sampleVelocity[i]);
i++;
}
}
@@ -461,24 +402,17 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
c.fieldIOobject("sampledPositions", IOobject::NO_READ),
np
);
-// vectorFieldIOField sampleVelocity
-// (
-// c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
-// np
-// );
label i = 0;
forAllConstIter(Cloud<streamLineParticle>, c, iter)
{
lifeTime[i] = iter().lifeTime_;
sampledPositions[i] = iter().sampledPositions_;
-// sampleVelocity[i] = iter().sampleVelocity_;
i++;
}
lifeTime.write();
sampledPositions.write();
-// sampleVelocity.write();
}
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.H b/src/functionObjects/field/streamLine/streamLineParticle.H
index bcb4fd68c2683b0f7963753b394e4b22840600f8..42282e1790d500ac3a4603f60df314b846e502d5 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/functionObjects/field/streamLine/streamLineParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -46,16 +46,12 @@ SourceFiles
namespace Foam
{
-class streamLineParticleCloud;
-
-
// Forward declaration of friend functions and operators
class streamLineParticle;
Ostream& operator<<(Ostream&, const streamLineParticle&);
-
/*---------------------------------------------------------------------------*\
Class streamLineParticle Declaration
\*---------------------------------------------------------------------------*/
@@ -64,32 +60,38 @@ class streamLineParticle
:
public particle
{
-
public:
- //- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<Cloud<streamLineParticle>>
{
-
public:
+ // Public data
+
+ const PtrList<interpolation<scalar>>& vsInterp_;
+
+ const PtrList<interpolation<vector>>& vvInterp_;
+
+ const label UIndex_;
+
+ const bool trackForward_;
+
+ const label nSubCycle_;
- const PtrList<interpolation<scalar>>& vsInterp_;
- const PtrList<interpolation<vector>>& vvInterp_;
- const label UIndex_;
- const bool trackForward_;
- const label nSubCycle_;
- const scalar trackLength_;
+ const scalar trackLength_;
- DynamicList<vectorList>& allPositions_;
- List<DynamicList<scalarList>>& allScalars_;
- List<DynamicList<vectorList>>& allVectors_;
+ DynamicList<vectorList>& allPositions_;
+
+ List<DynamicList<scalarList>>& allScalars_;
+
+ List<DynamicList<vectorList>>& allVectors_;
// Constructors
+ //- Construct from components
trackingData
(
Cloud<streamLineParticle>& cloud,
@@ -99,7 +101,6 @@ public:
const bool trackForward,
const label nSubCycle,
const scalar trackLength,
-
DynamicList<List<point>>& allPositions,
List<DynamicList<scalarList>>& allScalars,
List<DynamicList<vectorList>>& allVectors
@@ -112,7 +113,6 @@ public:
trackForward_(trackForward),
nSubCycle_(nSubCycle),
trackLength_(trackLength),
-
allPositions_(allPositions),
allScalars_(allScalars),
allVectors_(allVectors)
@@ -139,22 +139,6 @@ private:
// 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.
vector interpolateFields
(
@@ -195,8 +179,7 @@ public:
return autoPtr<particle>(new streamLineParticle(*this));
}
- //- Factory class to read-construct particles used for
- // parallel transfer
+ //- Factory class to read-construct particles used for parallel transfer
class iNew
{
const polyMesh& mesh_;
@@ -223,8 +206,7 @@ public:
// Tracking
//- 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
// Executed before other patch-hitting functions
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index fadaf872bc1ab1f014cf4c011196de65ae952433..f16a25aa50e18366d84764f797fe7ae07f3117ff 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -176,16 +176,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
U += velocity;
- addNewParcel
- (
- p,
- U,
- Ei,
- celli,
- cellTetIs.face(),
- cellTetIs.tetPt(),
- typeId
- );
+ addNewParcel(p, celli, U, Ei, typeId);
}
}
}
@@ -465,27 +456,13 @@ template<class ParcelType>
void Foam::DSMCCloud<ParcelType>::addNewParcel
(
const vector& position,
+ const label celli,
const vector& U,
const scalar Ei,
- const label celli,
- const label tetFacei,
- const label tetPti,
const label typeId
)
{
- ParcelType* pPtr = new ParcelType
- (
- mesh_,
- position,
- U,
- Ei,
- celli,
- tetFacei,
- tetPti,
- typeId
- );
-
- this->addParticle(pPtr);
+ this->addParticle(new ParcelType(mesh_, position, celli, U, Ei, typeId));
}
@@ -1110,9 +1087,7 @@ void Foam::DSMCCloud<ParcelType>::dumpParticlePositions() const
template<class ParcelType>
void Foam::DSMCCloud<ParcelType>::autoMap(const mapPolyMesh& mapper)
{
- typedef typename ParcelType::trackingData tdType;
- tdType td(*this);
- Cloud<ParcelType>::template autoMap<tdType>(td, mapper);
+ Cloud<ParcelType>::autoMap(mapper);
// Update the cell occupancy field
cellOccupancy_.setSize(mesh_.nCells());
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
index d7ea2be9b4150d2e310b4bf3e67a1c5456a8900e..81bbc2ff649131997f1086be1ff1d1cc37e6f72f 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -44,6 +44,7 @@ SourceFiles
#include "fvMesh.H"
#include "volFields.H"
#include "scalarIOField.H"
+#include "barycentric.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -451,11 +452,9 @@ public:
void addNewParcel
(
const vector& position,
+ const label celli,
const vector& U,
const scalar Ei,
- const label celli,
- const label tetFacei,
- const label tetPti,
const label typeId
);
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
index 554f3b52bad479f1a6bf6ea772e76b74451d62b0..3df37ae16068f3d81b2216a8c8b60dbd27ea1298 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -41,19 +41,16 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
const polyMesh& mesh = td.cloud().pMesh();
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
// constrained, but the actual U_ of the parcel must not be
// altered or used, as it is altered by patch interactions an
// needs to retain its 3D value for collision purposes.
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
- meshTools::constrainToMeshCentre(mesh, p.position());
+ p.constrainToMeshCentre();
Utracking = U_;
@@ -61,16 +58,10 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
// reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
- // Set the Lagrangian time-step
- scalar dt = min(dtMax, tEnd);
-
- dt *= p.trackToFace(p.position() + dt*Utracking, td);
-
- tEnd -= dt;
-
- p.stepFraction() = 1.0 - tEnd/trackTime;
+ const scalar f = 1 - p.stepFraction();
+ p.trackToFace(f*trackTime*Utracking, f, td);
- if (p.onBoundary() && td.keepParticle)
+ if (p.onBoundaryFace() && td.keepParticle)
{
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
index 453cb74315f3fb964dc7cac2b8a77e8e1c36e494..b3525cd7db4a6d1c15cc71f3bfb82608b654758e 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -176,12 +176,24 @@ public:
inline DSMCParcel
(
const polyMesh& mesh,
- const vector& position,
- const vector& U,
- const scalar Ei,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
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
);
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
index 5fb82c5cdad57bac557e87d555275b67f5d70f39..2bff944e14878ef58eb0f524df6024e7ab6196ec 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -55,16 +55,34 @@ template<class ParcelType>
inline Foam::DSMCParcel<ParcelType>::DSMCParcel
(
const polyMesh& mesh,
- const vector& position,
- const vector& U,
- const scalar Ei,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
+ const vector& U,
+ const scalar Ei,
+ const label typeId
+)
+:
+ 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, tetFacei, tetPti),
+ ParcelType(mesh, position, celli),
U_(U),
Ei_(Ei),
typeId_(typeId)
diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 4ed53b5cf68a9d901ad8845282c8559e9614527d..3b85d705876ba385ed5f54fb2a0d3eee13afd11b 100644
--- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -403,16 +403,7 @@ void Foam::FreeStream<CloudType>::inflow()
cloud.constProps(typeId).internalDegreesOfFreedom()
);
- cloud.addNewParcel
- (
- p,
- U,
- Ei,
- celli,
- globalFaceIndex,
- faceTetIs.tetPt(),
- typeId
- );
+ cloud.addNewParcel(p, celli, U, Ei, typeId);
particlesInserted++;
}
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 9b66f53706aabc512eba6538d2b50e7fdad4c8cb..f0f7da835c4db340e8b996d1aa55e519021fa561 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -105,7 +105,8 @@ Foam::Cloud<ParticleType>::Cloud
polyMesh_(pMesh),
labels_(),
nTrackingRescues_(),
- cellWallFacesPtr_()
+ cellWallFacesPtr_(),
+ globalPositionsPtr_()
{
checkPatches();
@@ -234,6 +235,8 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
// Allocate transfer buffers
PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+ // Clear the global positions as there are about to change
+ globalPositionsPtr_.clear();
// While there are particles to transfer
while (true)
@@ -393,21 +396,16 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
template<class ParticleType>
-template<class TrackData>
-void Foam::Cloud<ParticleType>::autoMap
-(
- TrackData& td,
- const mapPolyMesh& mapper
-)
+void Foam::Cloud<ParticleType>::autoMap(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
// polyMesh_.clearCellTree();
cellWallFacesPtr_.clear();
@@ -417,51 +415,13 @@ void Foam::Cloud<ParticleType>::autoMap
// there is a comms mismatch.
polyMesh_.tetBasePtIs();
+ const vectorField& positions = globalPositionsPtr_();
- forAllIter(typename Cloud<ParticleType>, *this, pIter)
+ label i = 0;
+ forAllIter(typename Cloud<ParticleType>, *this, iter)
{
- ParticleType& p = pIter();
-
- 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);
- }
+ iter().autoMap(positions[i], mapper);
+ ++ i;
}
}
@@ -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 * * * * * * * * * * * //
#include "CloudIO.C"
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index 2ad874691b7a492482d8192eed4a3ec89d64fba9..27be0c3a69745c14b01abeeaec8d06c5d67b7e50 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -87,6 +87,9 @@ class Cloud
//- Does the cell have wall faces
mutable autoPtr<PackedBoolList> cellWallFacesPtr_;
+ //- Temporary storage for the global particle positions
+ mutable autoPtr<vectorField> globalPositionsPtr_;
+
// Private Member Functions
@@ -162,7 +165,7 @@ public:
return IDLList<ParticleType>::size();
};
- DynamicList<label>& labels()
+ DynamicList<label>& labels() const
{
return labels_;
}
@@ -256,8 +259,7 @@ public:
//- Remap the cells of particles corresponding to the
// mesh topology change
- template<class TrackData>
- void autoMap(TrackData& td, const mapPolyMesh&);
+ void autoMap(const mapPolyMesh&);
// Read
@@ -304,6 +306,10 @@ public:
//- Write positions to \<cloudName\>_positions.obj file
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
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index 83369d1a4f0c28cca547af0324f4ded8a5417275..7e5131fd9b956db9d566a60bb1b7bc4affb5fd4d 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -132,13 +132,6 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
// them, otherwise, if some processors have no particles then
// there is a comms mismatch.
polyMesh_.tetBasePtIs();
-
- forAllIter(typename Cloud<ParticleType>, *this, pIter)
- {
- ParticleType& p = pIter();
-
- p.initCellFacePt();
- }
}
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 4223cd67789b672f81c3232fd5e466efe041fd02..b223232bd3ef7f8f34afc278407f14eca6f2dcad 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -965,14 +965,7 @@ void Foam::InteractionLists<ParticleType>::prepareParticleToBeReferred
globalTransforms.transformIndex(ciat)
);
- particle->position() = transform.invTransformPosition(particle->position());
-
- particle->transformProperties(-transform.t());
-
- if (transform.hasR())
- {
- particle->transformProperties(transform.R().T());
- }
+ particle->prepareForInteractionListReferral(transform);
}
@@ -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();
wallFaceMap().receive(pBufs, referredWallData_);
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticle.H b/src/lagrangian/basic/indexedParticle/indexedParticle.H
index 75ed0783e4284b142dd1e049935ebdaf5a73461e..c619cc6694a296644591b078379a302d4da9be49 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticle.H
+++ b/src/lagrangian/basic/indexedParticle/indexedParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -65,27 +65,14 @@ public:
indexedParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
const label index = 0
)
:
- particle(mesh, position, 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),
+ particle(mesh, coordinates, celli, tetFacei, tetPti),
index_(index)
{}
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 88de1e7dda789e9adc773d87f9fa9b17bc7e6af2..ea1a1a07a32554d7357c40bd69d32ee8d74d7a53 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,20 +25,549 @@ License
#include "particle.H"
#include "transform.H"
+#include "treeDataCell.H"
+#include "cubicEqn.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::label Foam::particle::particleCount_ = 0;
-const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
+namespace Foam
+{
+ defineTypeNameAndDebug(particle, 0);
+}
-const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
-const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-namespace Foam
+void Foam::particle::tetFaceIndices
+(
+ label& baseI,
+ label& vertex1I,
+ label& vertex2I
+) const
{
- defineTypeNameAndDebug(particle, 0);
+ const Foam::face& f = mesh_.faces()[tetFacei_];
+
+ baseI = max(0, mesh_.tetBasePtIs()[tetFacei_]);
+
+ vertex1I = (baseI + tetPti_) % f.size();
+
+ vertex2I = f.fcIndex(vertex1I);
+
+ if (mesh_.faceOwner()[tetFacei_] != celli_)
+ {
+ Swap(vertex1I, vertex2I);
+ }
+}
+
+
+void Foam::particle::tetMeshIndices
+(
+ label& basei,
+ label& vertex1i,
+ label& vertex2i
+) const
+{
+ const Foam::face& f = mesh_.faces()[tetFacei_];
+
+ tetFaceIndices(basei, vertex1i, vertex2i);
+
+ basei = f[basei];
+ vertex1i = f[vertex1i];
+ vertex2i = f[vertex2i];
+}
+
+
+void Foam::particle::tetGeometry
+(
+ vector& centre,
+ vector& base,
+ vector& vertex1,
+ vector& vertex2
+) const
+{
+ label basei, vertex1i, vertex2i;
+ tetMeshIndices(basei, vertex1i, vertex2i);
+
+ centre = mesh_.cellCentres()[celli_];
+ base = mesh_.points()[basei];
+ vertex1 = mesh_.points()[vertex1i];
+ vertex2 = mesh_.points()[vertex2i];
+}
+
+
+Foam::barycentricTensor Foam::particle::tetTransform() const
+{
+ vector centre, base, vertex1, vertex2;
+ tetGeometry(centre, base, vertex1, vertex2);
+
+ return barycentricTensor(centre, base, vertex1, vertex2);
+}
+
+
+void Foam::particle::tetReverseTransform
+(
+ vector& centre,
+ scalar& detA,
+ barycentricTensor& T
+) const
+{
+ barycentricTensor A = tetTransform();
+
+ const vector ab = A.b() - A.a();
+ const vector ac = A.c() - A.a();
+ const vector ad = A.d() - A.a();
+ const vector bc = A.c() - A.b();
+ const vector bd = A.d() - A.b();
+
+ centre = A.a();
+
+ detA = ab & (ac ^ ad);
+
+ T = barycentricTensor
+ (
+ bd ^ bc,
+ ac ^ ad,
+ ad ^ ab,
+ ab ^ ac
+ );
+}
+
+
+void Foam::particle::movingTetGeometry
+(
+ const scalar fraction,
+ Pair<vector>& centre,
+ Pair<vector>& base,
+ Pair<vector>& vertex1,
+ Pair<vector>& vertex2
+) const
+{
+ label basei, vertex1i, vertex2i;
+ tetMeshIndices(basei, vertex1i, vertex2i);
+
+ const pointField& ptsOld = mesh_.oldPoints();
+ const pointField& ptsNew = mesh_.points();
+
+ // !!! <-- We would be better off using mesh_.cellCentres() here. However,
+ // we need to put a mesh_.oldCellCentres() method in for this to work. The
+ // values obtained from the mesh and those obtained from the cell do not
+ // necessarily match. See mantis #1993.
+ const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
+ const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+
+ scalar f0 = stepFraction_, f1 = fraction;
+
+ // Old and new points and cell centres are not sub-cycled. If we are sub-
+ // cycling, then we have to account for the timestep change here by
+ // modifying the fractions that we take of the old and new geometry.
+ if (mesh_.time().subCycling())
+ {
+ const TimeState& tsNew = mesh_.time();
+ const TimeState& tsOld = mesh_.time().prevTimeState();
+
+ const scalar tFrac =
+ (
+ (tsNew.value() - tsNew.deltaTValue())
+ - (tsOld.value() - tsOld.deltaTValue())
+ )/tsOld.deltaTValue();
+
+ const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
+
+ f0 = tFrac + dtFrac*f0;
+ f1 = dtFrac*f1;
+ }
+
+ centre[0] = ccOld + f0*(ccNew - ccOld);
+ base[0] = ptsOld[basei] + f0*(ptsNew[basei] - ptsOld[basei]);
+ vertex1[0] = ptsOld[vertex1i] + f0*(ptsNew[vertex1i] - ptsOld[vertex1i]);
+ vertex2[0] = ptsOld[vertex2i] + f0*(ptsNew[vertex2i] - ptsOld[vertex2i]);
+
+ centre[1] = f1*(ccNew - ccOld);
+ base[1] = f1*(ptsNew[basei] - ptsOld[basei]);
+ vertex1[1] = f1*(ptsNew[vertex1i] - ptsOld[vertex1i]);
+ vertex2[1] = f1*(ptsNew[vertex2i] - ptsOld[vertex2i]);
+}
+
+
+Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
+(
+ const scalar fraction
+) const
+{
+ Pair<vector> centre, base, vertex1, vertex2;
+ movingTetGeometry(fraction, centre, base, vertex1, vertex2);
+
+ return
+ Pair<barycentricTensor>
+ (
+ barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
+ barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
+ );
+}
+
+
+void Foam::particle::movingTetReverseTransform
+(
+ const scalar fraction,
+ Pair<vector>& centre,
+ FixedList<scalar, 4>& detA,
+ FixedList<barycentricTensor, 3>& T
+) const
+{
+ Pair<barycentricTensor> A = movingTetTransform(fraction);
+
+ const Pair<vector> ab(A[0].b() - A[0].a(), A[1].b() - A[1].a());
+ const Pair<vector> ac(A[0].c() - A[0].a(), A[1].c() - A[1].a());
+ const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
+ const Pair<vector> bc(A[0].c() - A[0].b(), A[1].c() - A[1].b());
+ const Pair<vector> bd(A[0].d() - A[0].b(), A[1].d() - A[1].b());
+
+ centre[0] = A[0].a();
+ centre[1] = A[1].a();
+
+ detA[0] = ab[0] & (ac[0] ^ ad[0]);
+ detA[1] =
+ (ab[1] & (ac[0] ^ ad[0]))
+ + (ab[0] & (ac[1] ^ ad[0]))
+ + (ab[0] & (ac[0] ^ ad[1]));
+ detA[2] =
+ (ab[0] & (ac[1] ^ ad[1]))
+ + (ab[1] & (ac[0] ^ ad[1]))
+ + (ab[1] & (ac[1] ^ ad[0]));
+ detA[3] = ab[1] & (ac[1] ^ ad[1]);
+
+ T[0] = barycentricTensor
+ (
+ bd[0] ^ bc[0],
+ ac[0] ^ ad[0],
+ ad[0] ^ ab[0],
+ ab[0] ^ ac[0]
+ );
+ T[1] = barycentricTensor
+ (
+ (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
+ (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
+ (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
+ (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
+ );
+ T[2] = barycentricTensor
+ (
+ bd[1] ^ bc[1],
+ ac[1] ^ ad[1],
+ ad[1] ^ ab[1],
+ ab[1] ^ ac[1]
+ );
+}
+
+
+void Foam::particle::reflect()
+{
+ Swap(coordinates_.c(), coordinates_.d());
+}
+
+
+void Foam::particle::rotate(const bool reverse)
+{
+ if (!reverse)
+ {
+ scalar temp = coordinates_.b();
+ coordinates_.b() = coordinates_.c();
+ coordinates_.c() = coordinates_.d();
+ coordinates_.d() = temp;
+ }
+ else
+ {
+ scalar temp = coordinates_.d();
+ coordinates_.d() = coordinates_.c();
+ coordinates_.c() = coordinates_.b();
+ coordinates_.b() = temp;
+ }
+}
+
+
+void Foam::particle::changeTet(const label tetTriI)
+{
+ const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
+
+ const label firstTetPtI = 1;
+ const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
+
+ if (tetTriI == 1)
+ {
+ changeFace(tetTriI);
+ }
+ else if (tetTriI == 2)
+ {
+ if (isOwner)
+ {
+ if (tetPti_ == lastTetPtI)
+ {
+ changeFace(tetTriI);
+ }
+ else
+ {
+ reflect();
+ tetPti_ += 1;
+ }
+ }
+ else
+ {
+ if (tetPti_ == firstTetPtI)
+ {
+ changeFace(tetTriI);
+ }
+ else
+ {
+ reflect();
+ tetPti_ -= 1;
+ }
+ }
+ }
+ else if (tetTriI == 3)
+ {
+ if (isOwner)
+ {
+ if (tetPti_ == firstTetPtI)
+ {
+ changeFace(tetTriI);
+ }
+ else
+ {
+ reflect();
+ tetPti_ -= 1;
+ }
+ }
+ else
+ {
+ if (tetPti_ == lastTetPtI)
+ {
+ changeFace(tetTriI);
+ }
+ else
+ {
+ reflect();
+ tetPti_ += 1;
+ }
+ }
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Changing tet without changing cell should only happen when the "
+ << "track is on triangle 1, 2 or 3."
+ << exit(FatalError);
+ }
+}
+
+
+void Foam::particle::changeFace(const label tetTriI)
+{
+ // Get the tet topology
+ label basei, vertex1i, vertex2i;
+ tetMeshIndices(basei, vertex1i, vertex2i);
+
+ // Get the shared edge and the pre-rotation
+ edge sharedEdge;
+ if (tetTriI == 1)
+ {
+ sharedEdge = edge(vertex1i, vertex2i);
+ }
+ else if (tetTriI == 2)
+ {
+ sharedEdge = edge(vertex2i, basei);
+ }
+ else if (tetTriI == 3)
+ {
+ sharedEdge = edge(basei, vertex1i);
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Changing face without changing cell should only happen when the"
+ << " track is on triangle 1, 2 or 3."
+ << exit(FatalError);
+ }
+
+ // Find the face in the same cell that shares the edge, and the
+ // corresponding tetrahedra point
+ tetPti_ = -1;
+ forAll(mesh_.cells()[celli_], cellFaceI)
+ {
+ const label newFaceI = mesh_.cells()[celli_][cellFaceI];
+ const class face& newFace = mesh_.faces()[newFaceI];
+
+ // Exclude the current face
+ if (tetFacei_ == newFaceI)
+ {
+ continue;
+ }
+
+ // Loop over the edges, looking for the shared one
+ label edgeComp = 0;
+ label edgeI = 0;
+ for (; edgeI < newFace.size(); ++ edgeI)
+ {
+ edgeComp = edge::compare(sharedEdge, newFace.faceEdge(edgeI));
+
+ if (edgeComp != 0)
+ {
+ break;
+ }
+ }
+
+ // If the face does not contain the edge, then move on to the next face
+ if (edgeComp == 0)
+ {
+ continue;
+ }
+
+ // Make the edge index relative to the base point
+ const label newBaseI = max(0, mesh_.tetBasePtIs()[newFaceI]);
+ edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
+
+ // If the edge is next the base point (i.e., the index is 0 or n - 1),
+ // then we swap it for the adjacent edge. This new edge is opposite the
+ // base point, and defines the tet with the original edge in it.
+ edgeI = min(max(1, edgeI), newFace.size() - 2);
+
+ // Set the new face and tet point
+ tetFacei_ = newFaceI;
+ tetPti_ = edgeI;
+
+ // Exit the loop now that the the tet point has been found
+ break;
+ }
+
+ if (tetPti_ == -1)
+ {
+ FatalErrorInFunction
+ << "The search for an edge-connected face and tet-point failed."
+ << exit(FatalError);
+ }
+
+ // Pre-rotation puts the shared edge opposite the base of the tetrahedron
+ if (sharedEdge.otherVertex(vertex1i) == -1)
+ {
+ rotate(false);
+ }
+ else if (sharedEdge.otherVertex(vertex2i) == -1)
+ {
+ rotate(true);
+ }
+
+ // Update the new tet topology
+ tetMeshIndices(basei, vertex1i, vertex2i);
+
+ // Reflect to account for the change of triangle orientation on the new face
+ reflect();
+
+ // Post rotation puts the shared edge back in the correct location
+ if (sharedEdge.otherVertex(vertex1i) == -1)
+ {
+ rotate(true);
+ }
+ else if (sharedEdge.otherVertex(vertex2i) == -1)
+ {
+ rotate(false);
+ }
+}
+
+
+void Foam::particle::changeCell()
+{
+ // Set the cell to be the one on the other side of the face
+ const label ownerCellI = mesh_.faceOwner()[tetFacei_];
+ const bool isOwner = celli_ == ownerCellI;
+ celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
+
+ // Reflect to account for the change of triangle orientation in the new cell
+ reflect();
+}
+
+
+void Foam::particle::locate
+(
+ const vector& position,
+ const vector* direction,
+ const label celli,
+ const bool boundaryFail,
+ const string boundaryMsg
+)
+{
+ celli_ = celli;
+
+ // Find the cell, if it has not been given
+ if (celli_ < 0)
+ {
+ celli_ = mesh_.cellTree().findInside(position);
+ }
+ if (celli_ < 0)
+ {
+ FatalErrorInFunction
+ << "Cell not found for particle position " << position << "."
+ << exit(FatalError);
+ }
+
+ // Put the particle at the cell centre and in a random tet
+ coordinates_ = barycentric(1, 0, 0, 0);
+ tetFacei_ = mesh_.cells()[celli_][0];
+ tetPti_ = 1;
+ facei_ = -1;
+
+ // Track to the injection point
+ track(position - mesh_.cellCentres()[celli_], 0);
+ if (!onFace())
+ {
+ return;
+ }
+
+ // We hit a boundary ...
+ if (boundaryFail)
+ {
+ FatalErrorInFunction << boundaryMsg << exit(FatalError);
+ }
+ else
+ {
+ // Re-do the track, but this time do the bit tangential to the
+ // direction/patch first. This gets us as close as possible to the
+ // original path/position.
+
+ if (direction == nullptr)
+ {
+ const polyPatch& p = mesh_.boundaryMesh()[patch(facei_)];
+ direction = &p.faceNormals()[patchFace(patch(facei_), facei_)];
+ }
+
+ const vector n = *direction/mag(*direction);
+ const vector s = position - mesh_.cellCentres()[celli_];
+ const vector sN = (s & n)*n;
+ const vector sT = s - sN;
+
+ coordinates_ = barycentric(1, 0, 0, 0);
+ tetFacei_ = mesh_.cells()[celli_][0];
+ tetPti_ = 1;
+ facei_ = -1;
+
+ track(sT, 0);
+ track(sN, 0);
+
+ static label nWarnings = 0;
+ static const label maxNWarnings = 100;
+ if (nWarnings < maxNWarnings)
+ {
+ WarningInFunction << boundaryMsg << endl;
+ ++ nWarnings;
+ }
+ if (nWarnings == maxNWarnings)
+ {
+ WarningInFunction
+ << "Suppressing any further warnings about particles being "
+ << "located outside of the mesh." << endl;
+ ++ nWarnings;
+ }
+ }
}
@@ -47,19 +576,19 @@ namespace Foam
Foam::particle::particle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
mesh_(mesh),
- position_(position),
+ coordinates_(coordinates),
celli_(celli),
- facei_(-1),
- stepFraction_(0.0),
tetFacei_(tetFacei),
tetPti_(tetPti),
+ facei_(-1),
+ stepFraction_(0.0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{}
@@ -69,36 +598,39 @@ Foam::particle::particle
(
const polyMesh& mesh,
const vector& position,
- const label celli,
- bool doCellFacePt
+ const label celli
)
:
mesh_(mesh),
- position_(position),
+ coordinates_(- VGREAT, - VGREAT, - VGREAT, - VGREAT),
celli_(celli),
- facei_(-1),
- stepFraction_(0.0),
tetFacei_(-1),
tetPti_(-1),
+ facei_(-1),
+ stepFraction_(0.0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
- if (doCellFacePt)
- {
- initCellFacePt();
- }
+ locate
+ (
+ position,
+ nullptr,
+ celli,
+ false,
+ "Particle initialised with a location outside of the mesh."
+ );
}
Foam::particle::particle(const particle& p)
:
mesh_(p.mesh_),
- position_(p.position_),
+ coordinates_(p.coordinates_),
celli_(p.celli_),
- facei_(p.facei_),
- stepFraction_(p.stepFraction_),
tetFacei_(p.tetFacei_),
tetPti_(p.tetPti_),
+ facei_(p.facei_),
+ stepFraction_(p.stepFraction_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@@ -107,12 +639,12 @@ Foam::particle::particle(const particle& p)
Foam::particle::particle(const particle& p, const polyMesh& mesh)
:
mesh_(mesh),
- position_(p.position_),
+ coordinates_(p.coordinates_),
celli_(p.celli_),
- facei_(p.facei_),
- stepFraction_(p.stepFraction_),
tetFacei_(p.tetFacei_),
tetPti_(p.tetPti_),
+ facei_(p.facei_),
+ stepFraction_(p.stepFraction_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@@ -120,6 +652,319 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+Foam::scalar Foam::particle::track
+(
+ const vector& displacement,
+ const scalar fraction
+)
+{
+ scalar f = trackToFace(displacement, fraction);
+
+ while (onInternalFace())
+ {
+ changeCell();
+
+ f *= trackToFace(f*displacement, f*fraction);
+ }
+
+ return f;
+}
+
+
+Foam::scalar Foam::particle::trackToFace
+(
+ const vector& displacement,
+ const scalar fraction
+)
+{
+ scalar f = 1;
+
+ label tetTriI = onFace() ? 0 : -1;
+
+ facei_ = -1;
+
+ while (true)
+ {
+ if (mesh_.moving())
+ {
+ f *= trackToMovingTri(f*displacement, f*fraction, tetTriI);
+ }
+ else
+ {
+ f *= trackToStationaryTri(f*displacement, f*fraction, tetTriI);
+ }
+
+ if (tetTriI == -1)
+ {
+ // The track has completed within the current tet
+ return 0;
+ }
+ else if (tetTriI == 0)
+ {
+ // The track has hit a face, so set the current face and return
+ facei_ = tetFacei_;
+ return f;
+ }
+ else
+ {
+ // Move to the next tet and continue the track
+ changeTet(tetTriI);
+ }
+ }
+}
+
+
+Foam::scalar Foam::particle::trackToStationaryTri
+(
+ const vector& displacement,
+ const scalar fraction,
+ label& tetTriI
+)
+{
+ const vector x0 = position();
+ const vector x1 = displacement;
+ const barycentric y0 = coordinates_;
+
+ if (debug)
+ {
+ Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
+ }
+
+ // Get the tet geometry
+ vector centre;
+ scalar detA;
+ barycentricTensor T;
+ tetReverseTransform(centre, detA, T);
+
+ if (debug)
+ {
+ vector o, b, v1, v2;
+ tetGeometry(o, b, v1, v2);
+ Info<< "Tet points o=" << o << ", b=" << b
+ << ", v1=" << v1 << ", v2=" << v2 << endl
+ << "Tet determinant = " << detA << endl
+ << "Start local coordinates = " << y0 << endl;
+ }
+
+ // Calculate the local tracking displacement
+ barycentric Tx1(x1 & T);
+
+ if (debug)
+ {
+ Info<< "Local displacement = " << Tx1 << "/" << detA << endl;
+ }
+
+ // Calculate the hit fraction
+ label iH = -1;
+ scalar muH = detA == 0 ? VGREAT : mag(1/detA);
+ for (label i = 0; i < 4; ++ i)
+ {
+ if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
+ {
+ scalar mu = - y0[i]/Tx1[i];
+
+ if (debug)
+ {
+ Info<< "Hit on tet face " << i << " at local coordinate "
+ << y0 + mu*Tx1 << ", " << mu*detA*100 << "\% of the "
+ << "way along the track" << endl;
+ }
+
+ if (0 <= mu && mu < muH)
+ {
+ iH = i;
+ muH = mu;
+ }
+ }
+ }
+
+ // Set the new coordinates
+ barycentric yH = y0 + muH*Tx1;
+
+ // Clamp to zero any negative coordinates generated by round-off error
+ for (label i = 0; i < 4; ++ i)
+ {
+ yH.replace(i, i == iH ? 0 : max(0, yH[i]));
+ }
+
+ // Re-normalise if within the tet
+ if (iH == -1)
+ {
+ yH /= cmptSum(yH);
+ }
+
+ // Set the new position and hit index
+ coordinates_ = yH;
+ tetTriI = iH;
+
+ if (debug)
+ {
+ if (iH != -1)
+ {
+ Info<< "Track hit tet face " << iH << " first" << endl;
+ }
+ else
+ {
+ Info<< "Track hit no tet faces" << endl;
+ }
+ Info<< "End local coordinates = " << yH << endl
+ << "End global coordinates = " << position() << endl
+ << "Tracking displacement = " << position() - x0 << endl
+ << muH*detA*100 << "\% of the step from " << stepFraction_ << " to "
+ << stepFraction_ + fraction << " completed" << endl << endl;
+ }
+
+ // Set the proportion of the track that has been completed
+ stepFraction_ += fraction*muH*detA;
+
+ return 1 - muH*detA;
+}
+
+
+Foam::scalar Foam::particle::trackToMovingTri
+(
+ const vector& displacement,
+ const scalar fraction,
+ label& tetTriI
+)
+{
+ const vector x0 = position();
+ const vector x1 = displacement;
+ const barycentric y0 = coordinates_;
+
+ if (debug)
+ {
+ Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
+ }
+
+ // Get the tet geometry
+ Pair<vector> centre;
+ FixedList<scalar, 4> detA;
+ FixedList<barycentricTensor, 3> T;
+ movingTetReverseTransform(fraction, centre, detA, T);
+
+ const vector x0Rel = x0 - centre[0];
+ const vector x1Rel = x1 - centre[1];
+
+ // Form the determinant and hit equations
+ cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
+ const barycentric yC(1, 0, 0, 0);
+ const barycentric hitEqnA = ((x1Rel & T[2]) + detA[3]*yC)*sqr(detA[0]);
+ const barycentric hitEqnB =
+ ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
+ const barycentric hitEqnC = (x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC;
+ const barycentric hitEqnD = y0;
+ FixedList<cubicEqn, 4> hitEqn;
+ forAll(hitEqn, i)
+ {
+ hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
+ }
+
+ // Calculate the hit fraction
+ label iH = -1;
+ scalar muH = detA[0] == 0 ? VGREAT : mag(1/detA[0]);
+ for (label i = 0; i < 4; ++ i)
+ {
+ const Roots<3> mu = hitEqn[i].roots();
+
+ for (label j = 0; j < 3; ++ j)
+ {
+ if (mu.type(j) == roots::real && hitEqn[i].derivative(mu[j]) < 0)
+ {
+ if (0 <= mu[j] && mu[j] < muH)
+ {
+ iH = i;
+ muH = mu[j];
+ }
+ }
+ }
+ }
+
+ // Set the new coordinates
+ barycentric yH
+ (
+ hitEqn[0].value(muH),
+ hitEqn[1].value(muH),
+ hitEqn[2].value(muH),
+ hitEqn[3].value(muH)
+ );
+ // !!! <-- This fails if the tet collapses onto the particle, as detA tends
+ // to zero at the hit. In this instance, we can differentiate the hit and
+ // detA polynomials to find a limiting location, but this will not be on a
+ // triangle. We will then need to do a second track through the degenerate
+ // tet to find the final hit position. This second track is over zero
+ // distance and therefore can be of the static mesh type. This has not yet
+ // been implemented.
+ const scalar detAH = detAEqn.value(muH);
+ if (!std::isnormal(detAH))
+ {
+ FatalErrorInFunction
+ << "A moving tet collapsed onto a particle. This is not supported. "
+ << "The mesh is too poor, or the motion too severe, for particle "
+ << "tracking to function." << exit(FatalError);
+ }
+ yH /= detAH;
+
+ // Clamp to zero any negative coordinates generated by round-off error
+ for (label i = 0; i < 4; ++ i)
+ {
+ yH.replace(i, i == iH ? 0 : max(0, yH[i]));
+ }
+
+ // Re-normalise if within the tet
+ if (iH == -1)
+ {
+ yH /= cmptSum(yH);
+ }
+
+ // Set the new position and hit index
+ coordinates_ = yH;
+ tetTriI = iH;
+
+ // Set the proportion of the track that has been completed
+ stepFraction_ += fraction*muH*detA[0];
+
+ if (debug)
+ {
+ if (iH != -1)
+ {
+ Info<< "Track hit tet face " << iH << " first" << endl;
+ }
+ else
+ {
+ Info<< "Track hit no tet faces" << endl;
+ }
+ Info<< "End local coordinates = " << yH << endl
+ << "End global coordinates = " << position() << endl;
+ }
+
+ return 1 - muH*detA[0];
+}
+
+
+void Foam::particle::constrainToMeshCentre()
+{
+ const Vector<label>& dirs = mesh_.geometricD();
+
+ bool isConstrained = false;
+ forAll(dirs, dirI)
+ {
+ if (dirs[dirI] == -1)
+ {
+ isConstrained = true;
+ break;
+ }
+ }
+
+ if (isConstrained)
+ {
+ vector pos = position();
+ meshTools::constrainToMeshCentre(mesh_, pos);
+ track(pos - position(), 0);
+ }
+}
+
+
void Foam::particle::transformProperties(const tensor&)
{}
@@ -134,6 +979,107 @@ Foam::scalar Foam::particle::wallImpactDistance(const vector&) const
}
+void Foam::particle::prepareForInteractionListReferral
+(
+ const vectorTensorTransform& transform
+)
+{
+ // Get the transformed position
+ const vector pos = transform.invTransformPosition(position());
+
+ // Break the topology
+ celli_ = -1;
+ tetFacei_ = -1;
+ tetPti_ = -1;
+ facei_ = -1;
+
+ // Store the position in the barycentric data
+ coordinates_ = barycentric(1 - cmptSum(pos), pos.x(), pos.y(), pos.z());
+
+ // Transform the properties
+ transformProperties(- transform.t());
+ if (transform.hasR())
+ {
+ transformProperties(transform.R().T());
+ }
+}
+
+
+void Foam::particle::correctAfterInteractionListReferral(const label celli)
+{
+ // Get the position from the barycentric data
+ const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
+
+ // Create some arbitrary topology for the supplied cell
+ celli_ = celli;
+ tetFacei_ = mesh_.cells()[celli_][0];
+ tetPti_ = 1;
+ facei_ = -1;
+
+ // Get the reverse transform and directly set the coordinates from the
+ // position. This isn't likely to be correct; the particle is probably not
+ // in this tet. It will, however, generate the correct vector when the
+ // position method is called. A referred particle should never be tracked,
+ // so this approximate topology is good enough. By using the nearby cell we
+ // minimise the error associated with the incorrect topology.
+ if (mesh_.moving())
+ {
+ NotImplemented;
+ }
+ else
+ {
+ vector centre;
+ scalar detA;
+ barycentricTensor T;
+ tetReverseTransform(centre, detA, T);
+ coordinates_ = barycentric(1, 0, 0, 0) + ((pos - centre) & T/detA);
+ }
+}
+
+
+Foam::label Foam::particle::procTetPt
+(
+ const polyMesh& procMesh,
+ const label procCell,
+ const label procTetFace
+) const
+{
+ // The tet point on the procMesh differs from the current tet point if the
+ // mesh and procMesh faces are of differing orientation. The change is the
+ // same as in particle::correctAfterParallelTransfer.
+
+ if
+ (
+ (mesh_.faceOwner()[tetFacei_] == celli_)
+ == (procMesh.faceOwner()[procTetFace] == procCell)
+ )
+ {
+ return tetPti_;
+ }
+ else
+ {
+ return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
+ }
+}
+
+
+void Foam::particle::autoMap
+(
+ const vector& position,
+ const mapPolyMesh& mapper
+)
+{
+ locate
+ (
+ position,
+ nullptr,
+ mapper.reverseCellMap()[celli_],
+ true,
+ "Particle mapped to a location outside of the mesh."
+ );
+}
+
+
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
bool Foam::operator==(const particle& pA, const particle& pB)
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index d94db0dd5614db6e9c7472669bf431d69d4317d5..ea678968ce23118d56c1c0e49bec1d0e956e33ca 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -33,6 +33,8 @@ Description
#define particle_H
#include "vector.H"
+#include "barycentric.H"
+#include "barycentricTensor.H"
#include "Cloud.H"
#include "IDLList.H"
#include "pointField.H"
@@ -42,6 +44,7 @@ Description
#include "FixedList.H"
#include "polyMeshTetDecomposition.H"
#include "particleMacros.H"
+#include "vectorTensorTransform.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -130,25 +133,19 @@ public:
};
-protected:
+private:
- // Protected data
+ // Private data
//- Reference to the polyMesh database
const polyMesh& mesh_;
- //- Position of particle
- vector position_;
+ //- Coordinates of particle
+ barycentric coordinates_;
//- Index of the cell it is in
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
// particle is in
label tetFacei_;
@@ -158,6 +155,12 @@ protected:
// point.
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
label origProc_;
@@ -165,74 +168,137 @@ protected:
label origId_;
+private:
+
// Private Member Functions
- //- Find the tet tri faces between position and tet centre
- 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;
+ // Tetrahedra functions
- //- Find the lambda value for the line to-from across the
- // given tri face, where p = from + lambda*(to - from)
- inline scalar tetLambda
- (
- const vector& from,
- const vector& to,
- const label triI,
- const vector& tetArea,
- const label tetPlaneBasePtI,
- const label celli,
- const label tetFacei,
- const label tetPti,
- const scalar tol
- ) const;
+ //- Get indices into the current face for the face-bound vertices of
+ // the current tet.
+ void tetFaceIndices
+ (
+ label& baseI,
+ label& vertex1I,
+ label& vertex2I
+ ) const;
- //- Find the lambda value for a moving tri face
- inline scalar movingTetLambda
- (
- const vector& from,
- const vector& to,
- const label triI,
- const vector& tetArea,
- const label tetPlaneBasePtI,
- const label celli,
- const label tetFacei,
- const label tetPti,
- const scalar tol
- ) const;
+ //- Get indices into the mesh points for the face-bound vertices of
+ // the current tet.
+ void tetMeshIndices
+ (
+ label& baseI,
+ label& vertex1I,
+ label& vertex2I
+ ) const;
- //- Modify the tet owner data by crossing triI
- inline void tetNeighbour(label triI);
+ //- Get the vertices of the current tet
+ void tetGeometry
+ (
+ vector& centre,
+ vector& base,
+ vector& vertex1,
+ vector& vertex2
+ ) const;
- //- Cross the from the given face across the given edge of the
- // given cell to find the resulting face and tetPti
- inline void crossEdgeConnectedFace
- (
- const label& celli,
- label& tetFacei,
- label& tetPti,
- const edge& e
- );
+ //- Get the transformation associated with the current tet. This
+ // will convert a barycentric position within the tet to a
+ // cartesian position in the global coordinate system. The
+ // conversion is x = A & y, where x is the cartesian position, y is
+ // the barycentric position and A is the transformation tensor.
+ barycentricTensor tetTransform() const;
+
+ //- Get the reverse transform associated with the current tet. The
+ // conversion is detA*y = (x - centre) & T. The variables x, y and
+ // centre have the same meaning as for the forward transform. T is
+ // the transposed inverse of the forward transform tensor, A,
+ // multiplied by its determinant, detA. This separation allows
+ // the barycentric tracking algorithm to function on inverted or
+ // degenerate tetrahedra.
+ void tetReverseTransform
+ (
+ vector& centre,
+ scalar& detA,
+ barycentricTensor& T
+ ) const;
- //- Hit wall faces in the current cell if the
- //- wallImpactDistance is non-zero. They may not be in
- //- Different tets to the current.
- template<class CloudType>
- inline void hitWallFaces
- (
- const CloudType& td,
- const vector& from,
- const vector& to,
- scalar& lambdaMin,
- tetIndices& closestTetIs
- );
+ //- 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
@@ -315,8 +381,8 @@ public:
//- String representation of properties
DefinePropertyList
(
- "(Px Py Pz) celli facei stepFraction "
- "tetFacei tetPti origProc origId"
+ "(Px Py Pz) celli tetFacei tetPti "
+ "facei stepFraction origProc origId"
);
//- String representation of property types
@@ -328,17 +394,6 @@ public:
//- Cumulative particle counter - used to provode unique ID
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
@@ -346,20 +401,19 @@ public:
particle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components, tetFacei_ and tetPti_ are not
- // supplied so they will be deduced by a search
+ //- Construct from a position and a cell, searching for the rest of the
+ // required topology
particle
(
const polyMesh& mesh,
const vector& position,
- const label celli,
- bool doCellFacePt = true
+ const label celli
);
//- Construct from Istream
@@ -412,27 +466,15 @@ public:
//- Return the mesh database
inline const polyMesh& mesh() const;
- //- Return current particle position
- inline const vector& position() const;
-
- //- Return current particle position
- inline vector& position();
-
- //- Return current cell particle is in
- inline label& cell();
+ //- Return current particle coordinates
+ inline const barycentric& coordinates() const;
//- Return current cell particle is in
inline label cell() const;
- //- Return current tet face particle is in
- inline label& tetFace();
-
//- Return current tet face particle is in
inline label tetFace() const;
- //- Return current tet face particle is in
- inline label& tetPt();
-
//- Return current tet face particle is in
inline label tetPt() const;
@@ -453,43 +495,32 @@ public:
// on oldPoints
inline vector oldNormal() const;
- //- Return current face particle is on otherwise -1
- inline label& face();
-
//- Return current face particle is on otherwise -1
inline label face() const;
- //- Return the impact model to be used, soft or hard (default).
- inline bool softImpact() const;
-
//- Return the particle current time
inline scalar currentTime() const;
// Check
- //- Check the stored cell value (setting if necessary) and
- // initialise the tetFace and tetPt values
- inline void initCellFacePt();
-
- //- Is the particle on the boundary/(or outside the domain)?
- inline bool onBoundary() const;
+ //- Is the particle on a face?
+ inline bool onFace() const;
- //- Is this global face an internal face?
- inline bool internalFace(const label facei) const;
+ //- Is the particle on an internal face?
+ inline bool onInternalFace() const;
- //- Is this global face a boundary face?
- inline bool boundaryFace(const label facei) const;
+ //- Is the particle on a boundary face?
+ inline bool onBoundaryFace() const;
- //- Which patch is particle on
+ //- Which patch is particle on // <-- !!!
inline label patch(const label facei) const;
- //- Which face of this patch is this particle on
- inline label patchFace
- (
- const label patchi,
- const label facei
- ) const;
+ //- Which face of this patch is this particle on // <-- !!!
+ inline label patchFace(const label patchi, const label facei) const;
+
+ //- Return current particle position
+ inline vector position() const;
//- Return the fraction of time-step completed
inline scalar& stepFraction();
@@ -512,28 +543,57 @@ public:
// Track
- //- Track particle to end of trajectory
- // or until it hits the boundary.
- // On entry 'stepFraction()' should be set to the fraction of the
- // time-step at which the tracking starts and on exit it contains
- // the fraction of the time-step completed.
- // Returns the boundary face index if the track stops at the
- // boundary, -1 otherwise.
- template<class TrackData>
- label track(const vector& endPosition, TrackData& td);
-
- //- 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.
+ //- Track along the displacement for a given fraction of the overall
+ // step. End when the track is complete, or when a boundary is hit.
+ // On exit, stepFraction_ will have been incremented to the current
+ // position, and facei_ will be set to the index of the boundary
+ // face that was hit, or -1 if the track completed within a cell.
+ // The proportion of the displacement still to be completed is
+ // returned.
+ scalar track
+ (
+ const vector& displacement,
+ const scalar fraction
+ );
+
+ //- As particle::track, but also stops on internal faces.
+ scalar trackToFace
+ (
+ const vector& displacement,
+ const scalar fraction
+ );
+
+ //- As particle::trackToFace, but also stops on tet triangles. On
+ // exit, tetTriI is set to the index of the tet triangle that was
+ // 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>
- scalar trackToFace(const vector& endPosition, TrackData& td);
+ void trackToFace
+ (
+ const vector& displacement,
+ const scalar fraction,
+ TrackData& td
+ );
- //- Return the index of the face to be used in the interpolation
- // routine
- inline label faceInterpolation() const;
+ //- Set the constrained components of the particle position to the
+ // mesh centre.
+ void constrainToMeshCentre();
// Transformations
@@ -564,6 +624,38 @@ public:
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
//- Read the fields associated with the owner cloud
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 5df1a22e92ee6ae41534e79015e7a19870556a89..074bd8a5b0f01328af7bd2e2c17cdeecdfa1ace8 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -26,541 +26,6 @@ License
#include "polyMesh.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 * * * * * * * * * * * * * //
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_;
-}
-
-
-inline Foam::vector& Foam::particle::position()
-{
- return position_;
+ return coordinates_;
}
@@ -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
{
return tetFacei_;
}
-inline Foam::label& Foam::particle::tetFace()
-{
- return tetFacei_;
-}
-
-
inline Foam::label Foam::particle::tetPt() const
{
return tetPti_;
}
-inline Foam::label& Foam::particle::tetPt()
-{
- return tetPti_;
-}
-
-
inline Foam::tetIndices Foam::particle::currentTetIndices() const
{
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_);
+}
+
+
+inline bool Foam::particle::onBoundaryFace() const
+{
+ return onFace() && !mesh_.isInternalFace(facei_);
+}
+
+
+inline Foam::vector Foam::particle::position() const
+{
+ if (mesh_.moving())
{
- mesh_.findCellFacePt
- (
- position_,
- celli_,
- tetFacei_,
- tetPti_
- );
-
- if (debug && celli_ == -1)
- {
- WarningInFunction
- << "cell, tetFace and tetPt search failure for position "
- << position_
- << endl;
- }
+ return movingTetTransform(0)[0] & coordinates_;
}
else
{
- mesh_.findTetFacePt(celli_, position_, tetFacei_, tetPti_);
-
- 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;
- }
- }
+ return tetTransform() & coordinates_;
}
}
-inline bool Foam::particle::onBoundary() const
-{
- return facei_ != -1 && facei_ >= mesh_.nInternalFaces();
-}
-
-
inline Foam::scalar& Foam::particle::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
{
- return
- 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);
+ return mesh_.time().value() + stepFraction_*mesh_.time().deltaTValue();
}
@@ -894,10 +191,4 @@ inline Foam::label Foam::particle::patchFace
}
-inline Foam::label Foam::particle::faceInterpolation() const
-{
- return facei_;
-}
-
-
// ************************************************************************* //
diff --git a/src/lagrangian/basic/particle/particleIO.C b/src/lagrangian/basic/particle/particleIO.C
index b9e40bc3d9228fb6da0a26b74f4787edf9ef41d3..e2842be40910e8b2c9092b418d0a1ab4d8cac8c4 100644
--- a/src/lagrangian/basic/particle/particleIO.C
+++ b/src/lagrangian/basic/particle/particleIO.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -33,12 +33,12 @@ Foam::string Foam::particle::propertyTypes_ = Foam::particle::propertyTypes();
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
(
- sizeof(particle) - offsetof(particle, position_)
+ sizeof(particle) - offsetof(particle, coordinates_)
);
@@ -47,38 +47,33 @@ const std::size_t Foam::particle::sizeofFields
Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
:
mesh_(mesh),
- position_(),
+ coordinates_(),
celli_(-1),
- facei_(-1),
- stepFraction_(0.0),
tetFacei_(-1),
tetPti_(-1),
+ facei_(-1),
+ stepFraction_(0.0),
origProc_(Pstream::myProcNo()),
origId_(-1)
{
if (is.format() == IOstream::ASCII)
{
- is >> position_ >> celli_;
+ is >> coordinates_ >> celli_ >> tetFacei_ >> tetPti_;
if (readFields)
{
- is >> facei_
- >> stepFraction_
- >> tetFacei_
- >> tetPti_
- >> origProc_
- >> origId_;
+ is >> facei_ >> stepFraction_ >> origProc_ >> origId_;
}
}
else
{
if (readFields)
{
- is.read(reinterpret_cast<char*>(&position_), sizeofFields);
+ is.read(reinterpret_cast<char*>(&coordinates_), sizeofFields_);
}
else
{
- is.read(reinterpret_cast<char*>(&position_), sizeofPosition_);
+ is.read(reinterpret_cast<char*>(&coordinates_), sizeofPosition_);
}
}
@@ -91,11 +86,14 @@ void Foam::particle::writePosition(Ostream& os) const
{
if (os.format() == IOstream::ASCII)
{
- os << position_ << token::SPACE << celli_;
+ os << coordinates_
+ << token::SPACE << celli_
+ << token::SPACE << tetFacei_
+ << token::SPACE << tetPti_;
}
else
{
- os.write(reinterpret_cast<const char*>(&position_), sizeofPosition_);
+ os.write(reinterpret_cast<const char*>(&coordinates_), sizeofPosition_);
}
// Check state of Ostream
@@ -107,12 +105,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
{
if (os.format() == IOstream::ASCII)
{
- os << p.position_
+ os << p.coordinates_
<< token::SPACE << p.celli_
- << token::SPACE << p.facei_
- << token::SPACE << p.stepFraction_
<< token::SPACE << p.tetFacei_
<< token::SPACE << p.tetPti_
+ << token::SPACE << p.facei_
+ << token::SPACE << p.stepFraction_
<< token::SPACE << p.origProc_
<< token::SPACE << p.origId_;
}
@@ -120,8 +118,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
{
os.write
(
- reinterpret_cast<const char*>(&p.position_),
- particle::sizeofFields
+ reinterpret_cast<const char*>(&p.coordinates_),
+ particle::sizeofFields_
);
}
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index 95f018665c965b74f259351c7b1731e7f05c5693..49e88257368f0f2649bb8d48aaa0f897cdec6f41 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -44,7 +44,7 @@ void Foam::particle::prepareForParallelTransfer
)
{
// 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 =
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())
{
const tensor& T =
@@ -85,41 +79,22 @@ void Foam::particle::correctAfterParallelTransfer
transformProperties(-s);
}
- tetFacei_ = facei_ + ppp.start();
-
- // Faces either side of a coupled patch have matched base indices,
- // tetPti is specified relative to the base point, already and
- // opposite circulation directions by design, so if the vertices
- // are:
- // source:
- // 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.
-
+ // Set the topology
+ celli_ = ppp.faceCells()[facei_];
+ facei_ += ppp.start();
+ tetFacei_ = facei_;
+ // Faces either side of a coupled patch are numbered in opposite directions
+ // as their normals both point away from their connected cells. The tet
+ // point therefore counts in the opposite direction from the base point.
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
- // Reset the face index for the next tracking operation
- if (stepFraction_ > (1.0 - SMALL))
- {
- stepFraction_ = 1.0;
- facei_ = -1;
- }
- else
- {
- facei_ += ppp.start();
- }
+ // Reflect to account for the change of triangle orientation in the new cell
+ reflect();
+
+ // Note that the position does not need transforming explicitly. The face-
+ // triangle on the receive patch is the transformation of the one on the
+ // send patch, so whilst the barycentric coordinates remain the same, the
+ // change of triangle implicitly transforms the position.
}
@@ -158,7 +133,6 @@ void Foam::particle::readFields(CloudType& c)
template<class CloudType>
void Foam::particle::writeFields(const CloudType& c)
{
- // Write the cloud position file
IOPosition<CloudType> ioP(c);
ioP.write();
@@ -169,7 +143,11 @@ void Foam::particle::writeFields(const CloudType& c)
c.fieldIOobject("origProcId", IOobject::NO_READ),
np
);
- IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
+ IOField<label> origId
+ (
+ c.fieldIOobject("origId", IOobject::NO_READ),
+ np
+ );
label i = 0;
forAllConstIter(typename CloudType, c, iter)
@@ -208,387 +186,55 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
template<class TrackData>
-Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
-{
- 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
+void Foam::particle::trackToFace
(
- const vector& endPosition,
+ const vector& displacement,
+ const scalar fraction,
TrackData& td
)
{
- typedef typename TrackData::cloudType cloudType;
- typedef typename cloudType::particleType particleType;
-
- cloudType& cloud = td.cloud();
-
-
- 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)
- {
- // 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();
- }
+ // Track
+ trackToFace(displacement, fraction);
- 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);
+ // If the track is complete, return
+ if (!onFace())
+ {
+ return;
+ }
+ // Hit face/patch processing
+ typedef typename TrackData::cloudType::particleType particleType;
particleType& p = static_cast<particleType&>(*this);
p.hitFace(td);
-
- if (internalFace(facei_))
+ if (onInternalFace())
{
- // Change tet ownership because a tri face has been crossed,
- // 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);
- }
+ changeCell();
}
else
{
label origFacei = facei_;
- label patchi = patch(facei_);
+ label patchi = mesh_.boundaryMesh().whichPatch(facei_);
+
+ // No action is taken for tetPti_ for tetFacei_ here. These are handled
+ // by the patch interaction call or later during processor transfer.
- // No action taken for tetPti_ for tetFacei_ here, handled by
- // patch interaction call or later during processor transfer.
+ const tetIndices faceHitTetIs =
+ polyMeshTetDecomposition::triangleTetIndices
+ (
+ mesh_,
+ tetFacei_,
+ celli_,
+ tetPti_
+ );
if
(
- !p.hitPatch
+ !p.hitPatch
(
mesh_.boundaryMesh()[patchi],
td,
patchi,
- trackFraction,
+ stepFraction(),
faceHitTetIs
)
)
@@ -596,7 +242,7 @@ Foam::scalar Foam::particle::trackToFace
// Did patch interaction model switch patches?
if (facei_ != origFacei)
{
- patchi = patch(facei_);
+ patchi = mesh_.boundaryMesh().whichPatch(facei_);
}
const polyPatch& patch = mesh_.boundaryMesh()[patchi];
@@ -635,7 +281,7 @@ Foam::scalar Foam::particle::trackToFace
(
static_cast<const cyclicAMIPolyPatch&>(patch),
td,
- endPosition - position_
+ displacement
);
}
else if (isA<processorPolyPatch>(patch))
@@ -658,290 +304,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;
- }
- }
- }
- }
- }
}
@@ -1017,22 +379,17 @@ void Foam::particle::hitCyclicPatch
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_];
-
- tetFacei_ = facei_;
-
- // See note in correctAfterParallelTransfer for tetPti_ addressing.
+ // See note in correctAfterParallelTransfer for tetPti addressing ...
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
- const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
- label patchFacei = receiveCpp.whichFace(facei_);
-
- // Now the particle is on the receiving side
-
- // Have patch transform the position
- receiveCpp.transformPosition(position_, patchFacei);
+ // Reflect to account for the change of triangle orientation in the new cell
+ reflect();
// Transform the properties
if (!receiveCpp.parallel())
@@ -1041,7 +398,7 @@ void Foam::particle::hitCyclicPatch
(
receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0]
- : receiveCpp.forwardT()[patchFacei]
+ : receiveCpp.forwardT()[receiveFacei]
);
transformProperties(T);
}
@@ -1051,7 +408,7 @@ void Foam::particle::hitCyclicPatch
(
(receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0]
- : receiveCpp.separation()[patchFacei]
+ : receiveCpp.separation()[receiveFacei]
);
transformProperties(-s);
}
@@ -1066,38 +423,41 @@ void Foam::particle::hitCyclicAMIPatch
const vector& direction
)
{
- const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
-
- // Patch face index on sending side
- label patchFacei = facei_ - cpp.start();
+ vector pos = position();
- // Patch face index on receiving side - also updates position
- patchFacei = cpp.pointFace(patchFacei, direction, position_);
+ const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
+ const label sendFacei = cpp.whichFace(facei_);
+ const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
- if (patchFacei < 0)
+ if (receiveFacei < 0)
{
- // If the patch face of the particle is not known assume that
- // the particle is lost and to be deleted
+ // If the patch face of the particle is not known assume that the
+ // particle is lost and mark it to be deleted.
td.keepParticle = false;
-
WarningInFunction
<< "Particle lost across " << cyclicAMIPolyPatch::typeName
<< " patches " << cpp.name() << " and " << receiveCpp.name()
- << " at position " << position_
- << endl;
+ << " at position " << pos << endl;
}
- // Convert face index into global numbering
- facei_ = patchFacei + receiveCpp.start();
+ // Set the topology
+ facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
- celli_ = mesh_.faceOwner()[facei_];
-
- tetFacei_ = facei_;
-
- // See note in correctAfterParallelTransfer for tetPti_ addressing.
- tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
+ // 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."
+ );
- // Now the particle is on the receiving side
+ // The particle must remain associated with a face for the tracking to
+ // register as incomplete
+ facei_ = tetFacei_;
// Transform the properties
if (!receiveCpp.parallel())
@@ -1106,7 +466,7 @@ void Foam::particle::hitCyclicAMIPatch
(
receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0]
- : receiveCpp.forwardT()[patchFacei]
+ : receiveCpp.forwardT()[receiveFacei]
);
transformProperties(T);
}
@@ -1116,7 +476,7 @@ void Foam::particle::hitCyclicAMIPatch
(
(receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0]
- : receiveCpp.separation()[patchFacei]
+ : receiveCpp.separation()[receiveFacei]
);
transformProperties(-s);
}
diff --git a/src/lagrangian/basic/passiveParticle/passiveParticle.H b/src/lagrangian/basic/passiveParticle/passiveParticle.H
index 0624039fdc532120d9ddb885cac76a78747accd8..35ef33c0ae02d4cec58bc8412ddd2836650f16f4 100644
--- a/src/lagrangian/basic/passiveParticle/passiveParticle.H
+++ b/src/lagrangian/basic/passiveParticle/passiveParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -61,26 +61,26 @@ public:
passiveParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
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
(
const polyMesh& mesh,
const vector& position,
- const label celli,
- bool doCellFacePt = true
+ const label celli
)
:
- particle(mesh, position, celli, doCellFacePt)
+ particle(mesh, position, celli)
{}
//- Construct from Istream
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 0dfb75094d7652b2dc07710b8b2896cfb7af4f2b..6b22c85869f8b27f8bc50bad55a9facf12e3dd0b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -318,7 +318,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
),
rndGen_(Pstream::myProcNo()),
cellOccupancyPtr_(),
- cellLengthScale_(cbrt(mesh_.V())),
+ cellLengthScale_(mag(cbrt(mesh_.V()))),
rho_(rho),
U_(U),
mu_(mu),
@@ -845,18 +845,14 @@ void Foam::KinematicCloud<CloudType>::updateMesh()
{
updateCellOccupancy();
injectors_.updateMesh();
- cellLengthScale_ = cbrt(mesh_.V());
+ cellLengthScale_ = mag(cbrt(mesh_.V()));
}
template<class CloudType>
void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{
- typedef typename particle::TrackingData<KinematicCloud<CloudType>> tdType;
-
- tdType td(*this);
-
- Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+ Cloud<parcelType>::autoMap(mapper);
updateMesh();
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 64476622555970dd166375d9fbb867bb89f51ff2..04862d198ecf29c397e3310bc5c7fb1f8fa022f6 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -333,11 +333,7 @@ void Foam::ReactingCloud<CloudType>::evolve()
template<class CloudType>
void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{
- typedef typename particle::TrackingData<ReactingCloud<CloudType>> tdType;
-
- tdType td(*this);
-
- Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+ Cloud<parcelType>::autoMap(mapper);
this->updateMesh();
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 21b03b25dcee0ca0e200654294b5ddf05388e022..3d9d76d60c6c40c74e045ab356777d1f3182fdec 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -259,12 +259,7 @@ void Foam::ReactingMultiphaseCloud<CloudType>::autoMap
const mapPolyMesh& mapper
)
{
- typedef typename particle::TrackingData<ReactingMultiphaseCloud<CloudType>>
- tdType;
-
- tdType td(*this);
-
- Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+ Cloud<parcelType>::autoMap(mapper);
this->updateMesh();
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index cde6a17b62579edf4ba730b2968716b164bb8598..45f683308bbc0f190ffab98e82e1ec44837739c2 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -480,11 +480,7 @@ void Foam::ThermoCloud<CloudType>::evolve()
template<class CloudType>
void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{
- typedef typename particle::TrackingData<ThermoCloud<CloudType>> tdType;
-
- tdType td(*this);
-
- Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+ Cloud<parcelType>::autoMap(mapper);
this->updateMesh();
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
index e1a8c7544e8f1811c448a8d7a5ca079bfaceb2c1..ed405f81c1bda91772e22443a3367dda8588116a 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -172,22 +172,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline CollidingParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline CollidingParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
index c40c4e5d54e6530cf35dfedb84078bfed2b66fb3..c817a75bde005e27b90fafa62369be31fc7040cf 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -63,13 +63,13 @@ template<class ParcelType>
inline Foam::CollidingParcel<ParcelType>::CollidingParcel
(
const polyMesh& owner,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(owner, position, celli, tetFacei, tetPti),
+ ParcelType(owner, coordinates, celli, tetFacei, tetPti),
f_(Zero),
angularMomentum_(Zero),
torque_(Zero),
@@ -82,6 +82,22 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
(
const polyMesh& owner,
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 tetFacei,
const label tetPti,
@@ -99,7 +115,7 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
ParcelType
(
owner,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index b593d4393bff82a69edb1d7a43b77512d065e29e..c2d11aee3cb1505ccb5f4ec5b5cc0789ff604a8b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -268,71 +268,50 @@ bool Foam::KinematicParcel<ParcelType>::move
const cloudSolution& solution = td.cloud().solution();
const scalarField& cellLengthScale = td.cloud().cellLengthScale();
- scalar tEnd = (1.0 - p.stepFraction())*trackTime;
- scalar dtMax = solution.deltaTMax(trackTime);
-
- bool tracking = true;
- label nTrackingStalled = 0;
-
- while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+ while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
{
// 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();
+ const label celli = p.cell();
+ const scalar sfrac = p.stepFraction();
- // Set the Lagrangian time-step
- scalar dt = min(dtMax, tEnd);
+ // Total displacement over the time-step
+ const vector s = trackTime*U_;
- // Cache the parcel current cell as this will change if a face is hit
- const label celli = p.cell();
+ // Cell length scale
+ const scalar l = cellLengthScale[p.cell()];
- const scalar magU = mag(U_);
- if (p.active() && tracking && (magU > ROOTVSMALL))
+ // 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;
- const scalar deltaLMax = solution.deltaLMax(cellLengthScale[celli]);
- const scalar dCorr = min(d, deltaLMax);
- dt *=
- dCorr/d
- *p.trackToFace(p.position() + dCorr*U_/magU, td);
+ // Track to the next face
+ p.trackToFace(f*s, f, td);
}
-
- tEnd -= dt;
-
- const scalar newStepFraction = 1.0 - tEnd/trackTime;
-
- if (tracking)
+ else
{
- if
- (
- mag(p.stepFraction() - newStepFraction)
- < particle::minStepFractionTol
- )
+ // Abandon the track, and move to the end of the sub-step. If the
+ // the mesh is moving, this will implicitly move the parcel.
+ if (mesh.moving())
{
- nTrackingStalled++;
-
- if (nTrackingStalled > maxTrackAttempts)
- {
- tracking = false;
- }
- }
- else
- {
- nTrackingStalled = 0;
+ WarningInFunction
+ << "Tracking was abandoned on a moving mesh. Parcels may "
+ << "move unphysically as a result." << endl;
}
+ p.stepFraction() += f;
}
- p.stepFraction() = newStepFraction;
-
- bool calcParcel = true;
- if (!tracking && solution.steadyState())
- {
- calcParcel = false;
- }
+ const scalar dt = (p.stepFraction() - sfrac)*trackTime;
// Avoid problems with extremely small timesteps
- if ((dt > ROOTVSMALL) && calcParcel)
+ if (!td.cloud().solution().steadyState() && dt > ROOTVSMALL)
{
// Update cell based properties
p.setCellValues(td, dt, celli);
@@ -345,7 +324,7 @@ bool Foam::KinematicParcel<ParcelType>::move
p.calc(td, dt, celli);
}
- if (p.onBoundary() && td.keepParticle)
+ if (p.onBoundaryFace() && td.keepParticle)
{
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index dd45972ca2707d46219882222e686ac2e485e837..06a1e879db9a4d2f47c254faa306f9ce5a8a0f7c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -340,22 +340,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline KinematicParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline KinematicParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
@@ -495,9 +504,6 @@ public:
// Helper functions
- //- Return the index of the face used in the interpolation routine
- inline label faceInterpolation() const;
-
//- Cell owner mass
inline scalar massCell(const label celli) const;
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index a0e60d35ad09ef55a94d6e3743a2f4857b8a704c..749c136883801317399630889ff818c1f5e78b00 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -73,13 +73,13 @@ template<class ParcelType>
inline Foam::KinematicParcel<ParcelType>::KinematicParcel
(
const polyMesh& owner,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(owner, position, celli, tetFacei, tetPti),
+ ParcelType(owner, coordinates, celli, tetFacei, tetPti),
active_(true),
typeId_(-1),
nParticle_(0),
@@ -101,6 +101,31 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
(
const polyMesh& owner,
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 tetFacei,
const label tetPti,
@@ -112,7 +137,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
const constantProperties& constProps
)
:
- ParcelType(owner, position, celli, tetFacei, tetPti),
+ ParcelType(owner, coordinates, celli, tetFacei, tetPti),
active_(true),
typeId_(typeId),
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>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell
(
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
index c932590a9ca29e7bb6c671dc098718f3dd7af3ff..9bae31a577ce636cd4dfd5ac869b3bbbb2de12de 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -190,22 +190,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline MPPICParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline MPPICParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
index 1af5110019e76e153ba6e9d04f2a6bfce14cf683..f4fd43e14b0154c7800e9f60f4341d5429bf3b53 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -29,13 +29,13 @@ template<class ParcelType>
inline Foam::MPPICParcel<ParcelType>::MPPICParcel
(
const polyMesh& owner,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(owner, position, celli, tetFacei, tetPti),
+ ParcelType(owner, coordinates, celli, tetFacei, tetPti),
UCorrect_(Zero)
{}
@@ -45,6 +45,19 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
(
const polyMesh& owner,
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 tetFacei,
const label tetPti,
@@ -60,7 +73,7 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
ParcelType
(
owner,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index c9525627d41bd0196b0973f62c28ea60157a990e..0aaaf25a2a8963dc4c132a3597f29e00753b7fb0 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -278,23 +278,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, position and topology
// Other properties initialised as null
inline ReactingMultiphaseParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
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
inline ReactingMultiphaseParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
index c64e29b262fc9737d9364982f5c2e50cfcee969b..f689e91a9e0144a5f35deafefece7b604e52bcb0 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -68,13 +68,13 @@ template<class ParcelType>
inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(mesh, position, celli, tetFacei, tetPti),
+ ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
YGas_(0),
YLiquid_(0),
YSolid_(0),
@@ -87,6 +87,22 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
(
const polyMesh& mesh,
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 tetFacei,
const label tetPti,
@@ -108,7 +124,7 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
ParcelType
(
mesh,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index c5463d92968fce8aebe8c323151a31db5763bb8b..6335ca9dbd6b85adcac1f189e9af02a690037f43 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -232,22 +232,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline ReactingParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline ReactingParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index a6359c400626e4cff76b588fb8ea56f67c7319b5..08699da9b876351f0dffd55f8f6b7de36e4fe9c3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -63,13 +63,13 @@ template<class ParcelType>
inline Foam::ReactingParcel<ParcelType>::ReactingParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(mesh, position, celli, tetFacei, tetPti),
+ ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
mass0_(0.0),
Y_(0),
pc_(0.0)
@@ -81,6 +81,21 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
(
const polyMesh& mesh,
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 tetFacei,
const label tetPti,
@@ -99,7 +114,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
ParcelType
(
mesh,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index c0a6777e06cfdfcc97b6197597970a75cc938578..4129e25835d7daebaafbd30af36aa4f809d530d0 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -288,22 +288,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline ThermoParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline ThermoParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index 9b2f365b4474c9036321e27b1d257ed5478fe03a..9c35a55a6f1edb22da92dc753d696f5ef9275eab 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -74,13 +74,13 @@ template<class ParcelType>
inline Foam::ThermoParcel<ParcelType>::ThermoParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(mesh, position, celli, tetFacei, tetPti),
+ ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
T_(0.0),
Cp_(0.0),
Tc_(0.0),
@@ -93,6 +93,22 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
(
const polyMesh& mesh,
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 tetFacei,
const label tetPti,
@@ -110,7 +126,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
ParcelType
(
mesh,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 9f496a63f1a7ca3bc18a6d70fbdd63bf4110b565..16568fcc9a893e9dbe37ee1e85eb91fc7db2a145 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -485,8 +485,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
meshTools::constrainToMeshCentre(mesh, pos);
// Create a new parcel
- parcelType* pPtr =
- new parcelType(mesh, pos, celli, tetFacei, tetPti);
+ parcelType* pPtr = new parcelType(mesh, pos, celli);
// Check/set new parcel thermo properties
cloud.setParcelThermoProperties(*pPtr, dt);
@@ -610,8 +609,7 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
meshTools::constrainToMeshCentre(mesh, pos);
// Create a new parcel
- parcelType* pPtr =
- new parcelType(mesh, pos, celli, tetFacei, tetPti);
+ parcelType* pPtr = new parcelType(mesh, pos, celli);
// Check/set new parcel thermo properties
cloud.setParcelThermoProperties(*pPtr, 0.0);
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
index cca62669394693fae46ad8dfce99cbf2368dc1fa..eb8fc4416fa91806d8660e9dc786554e4d17cd30 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -146,16 +146,6 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
{
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 =
max
(
@@ -166,14 +156,7 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
// Create a new parcel
parcelType* pPtr =
- new parcelType
- (
- this->owner().pMesh(),
- pos,
- celli,
- tetFacei,
- tetPti
- );
+ new parcelType(this->owner().pMesh(), pos, celli);
// Check/set new parcel thermo properties
td.cloud().setParcelThermoProperties(*pPtr, 0.0);
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 08237cb39879a4e0dfdd70fbb8258ec3cec345e8..0cb1d22242e71472e23d6f692ef7b2797f13bb8b 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -446,7 +446,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
}
// 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];
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index a5ba420f34a8f88f5025afefbd8754fc82347157..f565adbfac2ca6472bd8941a9893e8eb6be1aa4d 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -85,18 +85,10 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
{
// Leapfrog tracking part
- scalar tEnd = (1.0 - stepFraction())*trackTime;
- scalar dtMax = tEnd;
-
- while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+ while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
{
- // set the lagrangian time-step
- scalar dt = min(dtMax, tEnd);
-
- dt *= trackToFace(position() + dt*v_, td);
-
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/trackTime;
+ const scalar f = 1 - stepFraction();
+ trackToFace(f*trackTime*v_, f, td);
}
}
else if (td.part() == 2)
@@ -205,7 +197,7 @@ void Foam::molecule::transformProperties(const tensor& T)
rf_ = transform(T, rf_);
- sitePositions_ = position_ + (T & (sitePositions_ - position_));
+ sitePositions_ = position() + (T & (sitePositions_ - position()));
siteForces_ = T & siteForces_;
}
@@ -226,7 +218,7 @@ void Foam::molecule::transformProperties(const vector& separation)
void Foam::molecule::setSitePositions(const constantProperties& constProps)
{
- sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
+ sitePositions_ = position() + (Q_ & constProps.siteReferencePositions());
}
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 46ace15276bcd01a522b899217814233920c61d5..5fd058d473e0d6bf74451ff680b411cb92807bf6 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -238,7 +238,7 @@ public:
inline molecule
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
@@ -253,6 +253,24 @@ public:
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
molecule
(
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index 5c3058d0d4ebc0c373d1996f7bbbde44003f1dfa..b1c9a354b4d623c35f17aef8b5d4a2bd8e51c9aa 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -220,7 +220,7 @@ inline Foam::molecule::constantProperties::constantProperties
inline Foam::molecule::molecule
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
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),
v_(v),
a_(a),
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 5318c57867bb597184c8b640ddb9209181acb867..3ff6175c449668e8770abb19267b245dba30db20 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -746,17 +746,8 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition
);
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
+ const label cell =
+ mesh_.cellTree().findInside(globalPosition);
if (findIndex(zone, cell) != -1)
{
@@ -764,8 +755,6 @@ void Foam::moleculeCloud::initialiseMolecules
(
globalPosition,
cell,
- tetFace,
- tetPt,
id,
tethered,
temperature,
@@ -834,17 +823,11 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition
);
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
+ const label cell =
+ mesh_.cellTree().findInside
+ (
+ globalPosition
+ );
if (findIndex(zone, cell) != -1)
{
@@ -852,8 +835,6 @@ void Foam::moleculeCloud::initialiseMolecules
(
globalPosition,
cell,
- tetFace,
- tetPt,
id,
tethered,
temperature,
@@ -913,17 +894,11 @@ void Foam::moleculeCloud::initialiseMolecules
globalPosition
);
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
+ const label cell =
+ mesh_.cellTree().findInside
+ (
+ globalPosition
+ );
if (findIndex(zone, cell) != -1)
{
@@ -931,8 +906,6 @@ void Foam::moleculeCloud::initialiseMolecules
(
globalPosition,
cell,
- tetFace,
- tetPt,
id,
tethered,
temperature,
@@ -996,26 +969,12 @@ void Foam::moleculeCloud::createMolecule
(
const point& position,
label cell,
- label tetFace,
- label tetPt,
label id,
bool tethered,
scalar temperature,
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);
label special = 0;
@@ -1068,8 +1027,6 @@ void Foam::moleculeCloud::createMolecule
mesh_,
position,
cell,
- tetFace,
- tetPt,
Q,
v,
Zero,
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 35a109e940b8aab8e8510acc111321fd5c9fe467..4dda921e89cc8d12e5b06d531bc9701840855375 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -114,8 +114,6 @@ private:
(
const point& position,
label cell,
- label tetFace,
- label tetPt,
label id,
bool tethered,
scalar temperature,
diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C
index b94e99a5e15e0d87617c9c0884a369e95e2b4d30..ec189b4a141720b7d2bbb9802723ba9d5387d2c9 100644
--- a/src/lagrangian/solidParticle/solidParticle.C
+++ b/src/lagrangian/solidParticle/solidParticle.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -43,34 +43,27 @@ bool Foam::solidParticle::move
td.switchProcessor = false;
td.keepParticle = true;
- const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();
+ const polyBoundaryMesh& pbMesh = mesh().boundaryMesh();
- scalar tEnd = (1.0 - stepFraction())*trackTime;
- scalar dtMax = tEnd;
-
- while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+ while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
{
if (debug)
{
- Info<< "Time = " << mesh_.time().timeName()
+ Info<< "Time = " << mesh().time().timeName()
<< " trackTime = " << trackTime
- << " tEnd = " << tEnd
<< " steptFraction() = " << stepFraction() << endl;
}
- // set the lagrangian time-step
- scalar dt = min(dtMax, tEnd);
- // remember which cell the parcel is in
- // since this will change if a face is hit
- label celli = cell();
+ const label celli = cell();
+ const scalar sfrac = stepFraction();
- dt *= trackToFace(position() + dt*U_, td);
+ const scalar f = 1 - stepFraction();
+ trackToFace(f*trackTime*U_, f, td);
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/trackTime;
+ const scalar dt = (stepFraction() - sfrac)*trackTime;
- cellPointWeight cpw(mesh_, position(), celli, face());
+ cellPointWeight cpw(mesh(), position(), celli, face());
scalar rhoc = td.rhoInterp().interpolate(cpw);
vector Uc = td.UInterp().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);
- if (onBoundary() && td.keepParticle)
+ if (onBoundaryFace() && td.keepParticle)
{
if (isA<processorPolyPatch>(pbMesh[patch(face())]))
{
@@ -133,7 +126,7 @@ void Foam::solidParticle::hitWallPatch
const tetIndices& tetIs
)
{
- vector nw = tetIs.faceTri(mesh_).normal();
+ vector nw = tetIs.faceTri(mesh()).normal();
nw /= mag(nw);
scalar Un = U_ & nw;
diff --git a/src/lagrangian/solidParticle/solidParticle.H b/src/lagrangian/solidParticle/solidParticle.H
index f763abf33952dc8b49fb3ea777539938a2d273fe..bc03cd3b2e30fd593c8982315b6bcb2307c54397 100644
--- a/src/lagrangian/solidParticle/solidParticle.H
+++ b/src/lagrangian/solidParticle/solidParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -133,7 +133,7 @@ public:
inline solidParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/solidParticle/solidParticleI.H b/src/lagrangian/solidParticle/solidParticleI.H
index c1b6fdd4189f49f70abdb164686623af43edeb9b..a9a43591a3cac2822d93c0715c2f66aa22c64a34 100644
--- a/src/lagrangian/solidParticle/solidParticleI.H
+++ b/src/lagrangian/solidParticle/solidParticleI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -45,7 +45,7 @@ inline Foam::solidParticle::trackingData::trackingData
inline Foam::solidParticle::solidParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
@@ -53,7 +53,7 @@ inline Foam::solidParticle::solidParticle
const vector& U
)
:
- particle(mesh, position, celli, tetFacei, tetPti),
+ particle(mesh, coordinates, celli, tetFacei, tetPti),
d_(d),
U_(U)
{}
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
index 98bf7fcc22709242a1d82cf6dc52e19b707c428b..c5cd250bfef6eaf90e08095a5bbcd179ba0a04f5 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -219,22 +219,31 @@ public:
// Constructors
- //- Construct from owner, position, and cloud owner
+ //- Construct from mesh, coordinates and topology
// Other properties initialised as null
inline SprayParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
);
- //- Construct from components
+ //- 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
+ inline SprayParcel
+ (
+ const polyMesh& mesh,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti,
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
index d0f7391b8269e02ff027888542a57f60e42bd280..44a8714e2cd7e8ea1d49cf85ff0fb1ad102fb758 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -107,15 +107,15 @@ template<class ParcelType>
inline Foam::SprayParcel<ParcelType>::SprayParcel
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
- ParcelType(mesh, position, celli, tetFacei, tetPti),
+ ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
d0_(this->d()),
- position0_(position),
+ position0_(this->position()),
sigma_(0.0),
mu_(0.0),
liquidCore_(0.0),
@@ -135,6 +135,31 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
(
const polyMesh& mesh,
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 tetFacei,
const label tetPti,
@@ -162,7 +187,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
ParcelType
(
mesh,
- position,
+ coordinates,
celli,
tetFacei,
tetPti,
@@ -178,7 +203,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
constProps
),
d0_(d0),
- position0_(position),
+ position0_(this->position()),
sigma_(constProps.sigma0()),
mu_(constProps.mu0()),
liquidCore_(liquidCore),
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 5d69b69098d80ed4acec1db750b8a759d0ddef84..12de3f54e2f622870bf09490409bda4eef40e353 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -350,18 +350,7 @@ void Foam::meshRefinement::markFeatureCellLevel
{
const point& keepPoint = keepPoints[i];
- label cellI = -1;
- 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);
- }
+ const label celli = mesh_.cellTree().findInside(keepPoint);
if (cellI != -1)
{
@@ -404,10 +393,15 @@ void Foam::meshRefinement::markFeatureCellLevel
(
mesh_,
keepPoint,
+<<<<<<< HEAD
cellI,
tetFaceI,
tetPtI,
featureMesh.points()[pointI], // endpos
+=======
+ celli,
+ featureMesh.points()[pointi], // endpos
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
featureLevel, // level
featI, // featureMesh
pointI, // end point
@@ -449,10 +443,15 @@ void Foam::meshRefinement::markFeatureCellLevel
(
mesh_,
keepPoint,
+<<<<<<< HEAD
cellI,
tetFaceI,
tetPtI,
featureMesh.points()[pointI], // endpos
+=======
+ celli,
+ featureMesh.points()[pointi], // endpos
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
featureLevel, // level
featI, // featureMesh
pointI, // end point
@@ -545,8 +544,14 @@ void Foam::meshRefinement::markFeatureCellLevel
label otherPointI = e.otherVertex(pointI);
trackedParticle* tp(new trackedParticle(startTp));
+<<<<<<< HEAD
tp->end() = featureMesh.points()[otherPointI];
tp->j() = otherPointI;
+=======
+ tp->start() = tp->position();
+ tp->end() = featureMesh.points()[otherPointi];
+ tp->j() = otherPointi;
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
tp->k() = edgeI;
if (debug&meshRefinement::FEATURESEEDS)
@@ -605,8 +610,14 @@ void Foam::meshRefinement::markFeatureCellLevel
const edge& e = featureMesh.edges()[edgeI];
label otherPointI = e.otherVertex(pointI);
+<<<<<<< HEAD
tp.end() = featureMesh.points()[otherPointI];
tp.j() = otherPointI;
+=======
+ tp.start() = tp.position();
+ tp.end() = featureMesh.points()[otherPointi];
+ tp.j() = otherPointi;
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
tp.k() = edgeI;
keepParticle = true;
break;
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
index e4ba1c54092df9844d90129ea7efd1f55b57da5a..71f44569be1f88749230677b3241abaff50af72a 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -31,7 +31,7 @@ License
Foam::trackedParticle::trackedParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPtI,
@@ -42,7 +42,30 @@ Foam::trackedParticle::trackedParticle
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),
level_(level),
i_(i),
@@ -64,7 +87,7 @@ Foam::trackedParticle::trackedParticle
{
if (is.format() == IOstream::ASCII)
{
- is >> end_;
+ is >> start_ >> end_;
level_ = readLabel(is);
i_ = readLabel(is);
j_ = readLabel(is);
@@ -74,8 +97,8 @@ Foam::trackedParticle::trackedParticle
{
is.read
(
- reinterpret_cast<char*>(&end_),
- sizeof(end_) + sizeof(level_)
+ reinterpret_cast<char*>(&start_),
+ sizeof(start_) + sizeof(end_) + sizeof(level_)
+ sizeof(i_) + sizeof(j_) + sizeof(k_)
);
}
@@ -96,9 +119,8 @@ bool Foam::trackedParticle::move
td.switchProcessor = false;
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
// on a processor boundary. If the endpoint is on a processor face
@@ -111,18 +133,13 @@ bool Foam::trackedParticle::move
{
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.
td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_);
- dt *= trackToFace(end_, td);
-
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/trackTime;
+ const scalar f = 1 - stepFraction();
+ trackToFace(f*(end_ - start_), f, td);
}
}
@@ -248,6 +265,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
if (os.format() == IOstream::ASCII)
{
os << static_cast<const particle&>(p)
+ << token::SPACE << p.start_
<< token::SPACE << p.end_
<< token::SPACE << p.level_
<< token::SPACE << p.i_
@@ -259,8 +277,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
os << static_cast<const particle&>(p);
os.write
(
- reinterpret_cast<const char*>(&p.end_),
- sizeof(p.end_) + sizeof(p.level_)
+ reinterpret_cast<const char*>(&p.start_),
+ sizeof(p.start_) + sizeof(p.end_) + sizeof(p.level_)
+ sizeof(p.i_) + sizeof(p.j_) + sizeof(p.k_)
);
}
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
index faf551946146e204aacaf4011b8887e97b22a98b..02a30b3a23d502c7093a8cc2b55181fd2b6fab3b 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -64,6 +64,9 @@ class trackedParticle
{
// Private data
+ //- Start point to track from
+ point start_;
+
//- End point to track to
point end_;
@@ -120,7 +123,7 @@ public:
trackedParticle
(
const polyMesh& mesh,
- const vector& position,
+ const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPtI,
@@ -131,6 +134,20 @@ public:
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
trackedParticle
(
@@ -170,6 +187,12 @@ public:
// Member Functions
+ //- Point to track from
+ point& start()
+ {
+ return start_;
+ }
+
//- Point to track to
point& end()
{
diff --git a/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C b/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
index 2136819b6c4481b3380da7996e646c0b0b69bc15..c8e8be4c5e0be65801dac7f2393ba1d8c080c8bf 100644
--- a/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
+++ b/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -48,9 +48,7 @@ void Foam::reconstructLagrangianPositions
forAll(meshes, i)
{
const labelList& cellMap = cellProcAddressing[i];
-
- // faceProcAddressing not required currently.
- // const labelList& faceMap = faceProcAddressing[i];
+ const labelList& faceMap = faceProcAddressing[i];
Cloud<passiveParticle> lpi(meshes[i], cloudName, false);
@@ -58,18 +56,21 @@ void Foam::reconstructLagrangianPositions
{
const passiveParticle& ppi = iter();
- // // Inverting sign if necessary and subtracting 1 from
- // // faceProcAddressing
- // label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
+ const label mappedCell = cellMap[ppi.cell()];
+
+ // Inverting sign if necessary and subtracting 1 from
+ // faceProcAddressing
+ label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
lagrangianPositions.append
(
new passiveParticle
(
mesh,
- ppi.position(),
- cellMap[ppi.cell()],
- false
+ ppi.coordinates(),
+ mappedCell,
+ mappedTetFace,
+ ppi.procTetPt(mesh, mappedCell, mappedTetFace)
)
);
}
diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
index b7d217a35c5b94a7d0c6eb5f613df5002f64d204..429d1c084badceea59b7e9265f300c2d7bf2c135 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.C
+++ b/src/sampling/sampledSet/face/faceOnlySet.C
@@ -56,13 +56,15 @@ bool Foam::faceOnlySet::trackToBoundary
{
particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
- const point& trackPt = singleParticle.position();
+ point trackPt = singleParticle.position();
while(true)
{
point oldPoint = trackPt;
- singleParticle.trackToFace(end_, trackData);
+ singleParticle.trackToFace(end_ - start_, 0, trackData);
+
+ trackPt = singleParticle.position();
if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
{
@@ -78,7 +80,7 @@ bool Foam::faceOnlySet::trackToBoundary
// End reached
return false;
}
- else if (singleParticle.onBoundary())
+ else if (singleParticle.onBoundaryFace())
{
// Boundary reached
return true;
diff --git a/src/sampling/sampledSet/polyLine/polyLineSet.C b/src/sampling/sampledSet/polyLine/polyLineSet.C
index 816945ec5b5eb4cd8c237387d5335d3895af100c..358162e7a6ec2e0c74559c7d8c886768c48a944b 100644
--- a/src/sampling/sampledSet/polyLine/polyLineSet.C
+++ b/src/sampling/sampledSet/polyLine/polyLineSet.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -56,29 +56,16 @@ bool Foam::polyLineSet::trackToBoundary
{
particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
- // Alias
- const point& trackPt = singleParticle.position();
-
while (true)
{
// Local geometry info
const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
const scalar smallDist = mag(tol*offset);
- point oldPos = trackPt;
- label facei = -1;
- do
- {
- singleParticle.stepFraction() = 0;
- singleParticle.track(sampleCoords_[sampleI+1], trackData);
- }
- while
- (
- !singleParticle.onBoundary()
- && (mag(trackPt - oldPos) < smallDist)
- );
+ singleParticle.track(offset, 0);
+ const point trackPt = singleParticle.position();
- if (singleParticle.onBoundary())
+ if (singleParticle.onBoundaryFace())
{
//Info<< "trackToBoundary : reached boundary"
// << " trackPt:" << trackPt << endl;
@@ -94,7 +81,7 @@ bool Foam::polyLineSet::trackToBoundary
// << endl;
samplingPts.append(trackPt);
samplingCells.append(singleParticle.cell());
- samplingFaces.append(facei);
+ samplingFaces.append(singleParticle.face());
// trackPt is at sampleI+1
samplingCurveDist.append(1.0*(sampleI+1));
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 78694160f3577a9ec06f99141755a70e2c54949d..1ad1364f7e37d3ba116a28d012dfc12712e2bcb5 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -209,17 +209,25 @@ Foam::point Foam::sampledSet::pushIn
label 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)
{
newPosition = facePt;
- label trap(1.0/particle::trackingCorrectionTol + 1);
+ label trap(1.0/trackingCorrectionTol + 1);
label iterNo = 0;
do
{
- newPosition += particle::trackingCorrectionTol*(cC - facePt);
+ newPosition += trackingCorrectionTol*(cC - facePt);
mesh().findTetFacePt
(
diff --git a/src/sampling/sampledSet/uniform/uniformSet.C b/src/sampling/sampledSet/uniform/uniformSet.C
index c16502d214071ac7ad9f0574c7ec580a4234f398..c09403a71cd981a5bfe338e1dccb03ebb5aa0a5f 100644
--- a/src/sampling/sampledSet/uniform/uniformSet.C
+++ b/src/sampling/sampledSet/uniform/uniformSet.C
@@ -94,8 +94,7 @@ bool Foam::uniformSet::trackToBoundary
const vector smallVec = tol*offset;
const scalar smallDist = mag(smallVec);
- // Alias
- const point& trackPt = singleParticle.position();
+ point trackPt = singleParticle.position();
particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
@@ -153,32 +152,10 @@ bool Foam::uniformSet::trackToBoundary
<< " to:" << samplePt << endl;
}
- point oldPos = trackPt;
- label facei = -1;
- do
- {
- singleParticle.stepFraction() = 0;
- singleParticle.track(samplePt, trackData);
-
- if (debug)
- {
- 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)
- );
+ singleParticle.track(samplePt - trackPt, 0);
+ trackPt = singleParticle.position();
- if (singleParticle.onBoundary())
+ if (singleParticle.onBoundaryFace())
{
//Pout<< "trackToBoundary : reached boundary" << endl;
if (mag(trackPt - samplePt) < smallDist)
@@ -188,7 +165,7 @@ bool Foam::uniformSet::trackToBoundary
// Reached samplePt on boundary
samplingPts.append(trackPt);
samplingCells.append(singleParticle.cell());
- samplingFaces.append(facei);
+ samplingFaces.append(singleParticle.face());
samplingCurveDist.append(mag(trackPt - start_));
}
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
index a626f4d688372f87d2a73d1ece56ba06125112c6..c767564e02d759f5d20c0377e5a50be8cfa6feab 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
@@ -47,7 +47,6 @@ runTimeModifiable true;
functions
{
#include "streamLines"
- #include "wallBoundedStreamLines"
#include "cuttingPlane"
#include "forceCoeffs"
#include "ensightWrite"