diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files
index 1ea75cf3ffffd1db870f44c9fa74fdb58afd050f..ba19ff80a06e3d3c1aceb2b7ea0e0aee97ebc708 100644
--- a/src/postProcessing/functionObjects/field/Make/files
+++ b/src/postProcessing/functionObjects/field/Make/files
@@ -29,6 +29,11 @@ streamLine/streamLineParticle.C
streamLine/streamLineParticleCloud.C
streamLine/streamLineFunctionObject.C
+wallBoundedStreamLine/wallBoundedStreamLine.C
+wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
+wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
+
surfaceInterpolateFields/surfaceInterpolateFields.C
surfaceInterpolateFields/surfaceInterpolateFieldsFunctionObject.C
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..4a3f4fe29c407c9526ce20dc3bdd03eae0b0c9c2
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict
@@ -0,0 +1,144 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 2.0.0 |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application icoFoam;
+
+startFrom startTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 0.5;
+
+deltaT 0.005;
+
+writeControl timeStep;
+
+writeInterval 20;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 6;
+
+writeCompression off;
+
+timeFormat general;
+
+timePrecision 6;
+
+runTimeModifiable true;
+
+functions
+{
+ readFields
+ {
+ // Where to load it from (if not already in solver)
+ functionObjectLibs ("libfieldFunctionObjects.so");
+
+ type readFields;
+ fields (p U k);
+ }
+
+ near
+ {
+ // Where to load it from
+ functionObjectLibs ("libfieldFunctionObjects.so");
+
+ type nearWallFields;
+
+ // Output every
+ outputControl outputTime;
+ //outputInterval 1;
+
+ // Fields to be sampled. Per field original name and mapped field to
+ // create.
+ fields
+ (
+ (U UNear)
+ );
+
+ // Patches/groups to sample (regular expressions)
+ patches (motorBike);
+
+ // Distance to sample
+ distance 0.001;
+ }
+ streamLines
+ {
+ // Where to load it from (if not already in solver)
+ functionObjectLibs ("libfieldFunctionObjects.so");
+ type wallBoundedStreamLine;
+
+ // Output every
+ outputControl timeStep; //outputTime;
+ // outputInterval 10;
+
+ setFormat vtk; //gnuplot; //xmgr; //raw; //jplot;
+
+ // Velocity field to use for tracking.
+ UName UNear;
+
+ // Interpolation method. Default is cellPoint. See sampleDict.
+ //interpolationScheme pointMVC;
+
+ // Tracked forwards (+U) or backwards (-U)
+ trackForward true;
+
+ interpolationScheme cellPoint;
+
+ // Names of fields to sample. Should contain above velocity field!
+ fields (p U k UNear);
+
+ // Steps particles can travel before being removed
+ lifeTime 100;
+
+ // Cloud name to use
+ cloudName particleTracks;
+
+ // Seeding method. See the sampleSets in sampleDict.
+ seedSampleSet patchSeed; //cloud;//triSurfaceMeshPointSet;
+
+ uniformCoeffs
+ {
+ type uniform;
+ axis x; //distance;
+
+ start (0.0035 0.0999 0.0001);
+ end (0.0035 0.0999 0.0099);
+ nPoints 20;
+ }
+
+ cloudCoeffs
+ {
+ type cloud;
+ axis x; //distance;
+ points ((0.351516548679288 -0.0116085375585099 1.24));
+ }
+ patchSeedCoeffs
+ {
+ type patchSeed;//patchCloud; //cloud; //uniform;
+ patches (motorBike);
+ axis x; //distance;
+ maxPoints 20000;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
new file mode 100644
index 0000000000000000000000000000000000000000..4b6f6761f72e32790df85fa8ad0ee8dd98963496
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
@@ -0,0 +1,887 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Pstream.H"
+#include "functionObjectList.H"
+#include "wallBoundedStreamLine.H"
+#include "fvMesh.H"
+#include "wallBoundedStreamLineParticleCloud.H"
+#include "ReadFields.H"
+#include "meshSearch.H"
+#include "sampledSet.H"
+#include "globalIndex.H"
+#include "mapDistribute.H"
+#include "interpolationCellPoint.H"
+#include "PatchTools.H"
+#include "meshSearchMeshObject.H"
+#include "faceSet.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::wallBoundedStreamLine, 0);
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::wallBoundedStreamLine::wallPatch() const
+{
+ const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+ label nFaces = 0;
+
+ forAll(patches, patchI)
+ {
+ //if (!polyPatch::constraintType(patches[patchI].type()))
+ if (isA<wallPolyPatch>(patches[patchI]))
+ {
+ nFaces += patches[patchI].size();
+ }
+ }
+
+ labelList addressing(nFaces);
+
+ nFaces = 0;
+
+ forAll(patches, patchI)
+ {
+ //if (!polyPatch::constraintType(patches[patchI].type()))
+ if (isA<wallPolyPatch>(patches[patchI]))
+ {
+ const polyPatch& pp = patches[patchI];
+
+ forAll(pp, i)
+ {
+ addressing[nFaces++] = pp.start()+i;
+ }
+ }
+ }
+
+ return autoPtr<indirectPrimitivePatch>
+ (
+ new indirectPrimitivePatch
+ (
+ IndirectList<face>
+ (
+ mesh.faces(),
+ addressing
+ ),
+ mesh.points()
+ )
+ );
+}
+
+
+Foam::tetIndices Foam::wallBoundedStreamLine::findNearestTet
+(
+ const PackedBoolList& isWallPatch,
+ const point& seedPt,
+ const label cellI
+) const
+{
+ const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+ const cell& cFaces = mesh.cells()[cellI];
+
+ label minFaceI = -1;
+ label minTetPtI = -1;
+ scalar minDistSqr = sqr(GREAT);
+
+ forAll(cFaces, cFaceI)
+ {
+ label faceI = cFaces[cFaceI];
+
+ if (isWallPatch[faceI])
+ {
+ const face& f = mesh.faces()[faceI];
+ const label fp0 = mesh.tetBasePtIs()[faceI];
+ const point& basePoint = mesh.points()[f[fp0]];
+
+ label fp = f.fcIndex(fp0);
+ for (label i = 2; i < f.size(); i++)
+ {
+ const point& thisPoint = mesh.points()[f[fp]];
+ label nextFp = f.fcIndex(fp);
+ const point& nextPoint = mesh.points()[f[nextFp]];
+
+ const triPointRef tri(basePoint, thisPoint, nextPoint);
+
+ scalar d2 = magSqr(tri.centre() - seedPt);
+ if (d2 < minDistSqr)
+ {
+ minDistSqr = d2;
+ minFaceI = faceI;
+ minTetPtI = i-1;
+ }
+ fp = nextFp;
+ }
+ }
+ }
+
+ // Put particle in tet
+ return tetIndices
+ (
+ cellI,
+ minFaceI,
+ minTetPtI,
+ mesh
+ );
+}
+
+
+void Foam::wallBoundedStreamLine::track()
+{
+ const Time& runTime = obr_.time();
+ const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+
+ // Determine the 'wall' patches
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // These are the faces that need to be followed
+
+ autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
+ PackedBoolList isWallPatch(mesh.nFaces());
+ forAll(boundaryPatch().addressing(), i)
+ {
+ isWallPatch[boundaryPatch().addressing()[i]] = 1;
+ }
+
+
+
+ // Find nearest wall particle for the seedPoints
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ IDLList<wallBoundedStreamLineParticle> initialParticles;
+ wallBoundedStreamLineParticleCloud particles
+ (
+ mesh,
+ cloudName_,
+ initialParticles
+ );
+
+ {
+ // Get the seed points
+ // ~~~~~~~~~~~~~~~~~~~
+
+ const sampledSet& seedPoints = sampledSetPtr_();
+
+
+ forAll(seedPoints, i)
+ {
+ const point& seedPt = seedPoints[i];
+ label cellI = seedPoints.cells()[i];
+
+ tetIndices ids(findNearestTet(isWallPatch, seedPt, cellI));
+
+ if (ids.face() != -1 && isWallPatch[ids.face()])
+ {
+ //Pout<< "Seeding particle :" << nl
+ // << " seedPt:" << seedPt << nl
+ // << " face :" << ids.face() << nl
+ // << " at :" << mesh.faceCentres()[ids.face()] << nl
+ // << " cell :" << mesh.cellCentres()[ids.cell()] << nl
+ // << endl;
+
+ particles.addParticle
+ (
+ new wallBoundedStreamLineParticle
+ (
+ mesh,
+ ids.faceTri(mesh).centre(),
+ ids.cell(),
+ ids.face(), // tetFace
+ ids.tetPt(),
+ -1, // not on a mesh edge
+ -1, // not on a diagonal edge
+ lifeTime_ // lifetime
+ )
+ );
+ }
+ else
+ {
+ Pout<< type() << " : ignoring seed " << seedPt
+ << " since not in wall cell." << endl;
+ }
+ }
+ }
+
+ label nSeeds = returnReduce(particles.size(), sumOp<label>());
+
+ Info<< type() << " : seeded " << nSeeds << " particles." << endl;
+
+
+
+ // Read or lookup fields
+ PtrList<volScalarField> vsFlds;
+ PtrList<interpolation<scalar> > vsInterp;
+ PtrList<volVectorField> vvFlds;
+ PtrList<interpolation<vector> > vvInterp;
+
+ label UIndex = -1;
+
+ if (loadFromFiles_)
+ {
+ IOobjectList allObjects(mesh, runTime.timeName());
+
+ IOobjectList objects(2*fields_.size());
+ forAll(fields_, i)
+ {
+ objects.add(*allObjects[fields_[i]]);
+ }
+
+ ReadFields(mesh, objects, vsFlds);
+ vsInterp.setSize(vsFlds.size());
+ forAll(vsFlds, i)
+ {
+ vsInterp.set
+ (
+ i,
+ interpolation<scalar>::New
+ (
+ interpolationScheme_,
+ vsFlds[i]
+ )
+ );
+ }
+ ReadFields(mesh, objects, vvFlds);
+ vvInterp.setSize(vvFlds.size());
+ forAll(vvFlds, i)
+ {
+ vvInterp.set
+ (
+ i,
+ interpolation<vector>::New
+ (
+ interpolationScheme_,
+ vvFlds[i]
+ )
+ );
+ }
+ }
+ else
+ {
+ label nScalar = 0;
+ label nVector = 0;
+
+ forAll(fields_, i)
+ {
+ if (mesh.foundObject<volScalarField>(fields_[i]))
+ {
+ nScalar++;
+ }
+ else if (mesh.foundObject<volVectorField>(fields_[i]))
+ {
+ nVector++;
+ }
+ else
+ {
+ FatalErrorIn("wallBoundedStreamLine::execute()")
+ << "Cannot find field " << fields_[i] << endl
+ << "Valid scalar fields are:"
+ << mesh.names(volScalarField::typeName) << endl
+ << "Valid vector fields are:"
+ << mesh.names(volVectorField::typeName)
+ << exit(FatalError);
+ }
+ }
+ vsInterp.setSize(nScalar);
+ nScalar = 0;
+ vvInterp.setSize(nVector);
+ nVector = 0;
+
+ forAll(fields_, i)
+ {
+ if (mesh.foundObject<volScalarField>(fields_[i]))
+ {
+ const volScalarField& f = mesh.lookupObject<volScalarField>
+ (
+ fields_[i]
+ );
+ vsInterp.set
+ (
+ nScalar++,
+ interpolation<scalar>::New
+ (
+ interpolationScheme_,
+ f
+ )
+ );
+ }
+ else if (mesh.foundObject<volVectorField>(fields_[i]))
+ {
+ const volVectorField& f = mesh.lookupObject<volVectorField>
+ (
+ fields_[i]
+ );
+
+ if (f.name() == UName_)
+ {
+ UIndex = nVector;
+ }
+
+ vvInterp.set
+ (
+ nVector++,
+ interpolation<vector>::New
+ (
+ interpolationScheme_,
+ f
+ )
+ );
+ }
+ }
+ }
+
+ // Store the names
+ scalarNames_.setSize(vsInterp.size());
+ forAll(vsInterp, i)
+ {
+ scalarNames_[i] = vsInterp[i].psi().name();
+ }
+ vectorNames_.setSize(vvInterp.size());
+ forAll(vvInterp, i)
+ {
+ vectorNames_[i] = vvInterp[i].psi().name();
+ }
+
+ // Check that we know the index of U in the interpolators.
+
+ if (UIndex == -1)
+ {
+ FatalErrorIn("wallBoundedStreamLine::execute()")
+ << "Cannot find field to move particles with : " << UName_
+ << endl
+ << "This field has to be present in the sampled fields "
+ << fields_
+ << " and in the objectRegistry." << endl
+ << exit(FatalError);
+ }
+
+ // Sampled data
+ // ~~~~~~~~~~~~
+
+ // Size to maximum expected sizes.
+ allTracks_.clear();
+ allTracks_.setCapacity(nSeeds);
+ allScalars_.setSize(vsInterp.size());
+ forAll(allScalars_, i)
+ {
+ allScalars_[i].clear();
+ allScalars_[i].setCapacity(nSeeds);
+ }
+ allVectors_.setSize(vvInterp.size());
+ forAll(allVectors_, i)
+ {
+ allVectors_[i].clear();
+ allVectors_[i].setCapacity(nSeeds);
+ }
+
+
+ // additional particle info
+ wallBoundedStreamLineParticle::trackingData td
+ (
+ particles,
+ vsInterp,
+ vvInterp,
+ UIndex, // index of U in vvInterp
+ trackForward_, // track in +u direction?
+ isWallPatch, // which faces are to follow
+
+ allTracks_,
+ allScalars_,
+ allVectors_
+ );
+
+
+ // Set very large dt. Note: cannot use GREAT since 1/GREAT is SMALL
+ // which is a trigger value for the tracking...
+ const scalar trackTime = Foam::sqrt(GREAT);
+
+ // Track
+ particles.move(td, trackTime);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLine::wallBoundedStreamLine
+(
+ const word& name,
+ const objectRegistry& obr,
+ const dictionary& dict,
+ const bool loadFromFiles
+)
+:
+ dict_(dict),
+ name_(name),
+ obr_(obr),
+ loadFromFiles_(loadFromFiles),
+ active_(true)
+{
+ // Only active if a fvMesh is available
+ if (isA<fvMesh>(obr_))
+ {
+ read(dict_);
+ }
+ else
+ {
+ active_ = false;
+ WarningIn
+ (
+ "wallBoundedStreamLine::wallBoundedStreamLine\n"
+ "(\n"
+ "const word&,\n"
+ "const objectRegistry&,\n"
+ "const dictionary&,\n"
+ "const bool\n"
+ ")"
+ ) << "No fvMesh available, deactivating."
+ << nl << endl;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLine::~wallBoundedStreamLine()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::wallBoundedStreamLine::read(const dictionary& dict)
+{
+ if (active_)
+ {
+ //dict_ = dict;
+ dict.lookup("fields") >> fields_;
+ if (dict.found("UName"))
+ {
+ dict.lookup("UName") >> UName_;
+ }
+ else
+ {
+ UName_ = "U";
+ if (dict.found("U"))
+ {
+ IOWarningIn
+ (
+ "wallBoundedStreamLine::read(const dictionary&)",
+ dict
+ ) << "Using deprecated entry \"U\"."
+ << " Please use \"UName\" instead."
+ << endl;
+ dict.lookup("U") >> UName_;
+ }
+ }
+
+ if (findIndex(fields_, UName_) == -1)
+ {
+ FatalIOErrorIn
+ (
+ "wallBoundedStreamLine::read(const dictionary&)",
+ dict
+ ) << "Velocity field for tracking " << UName_
+ << " should be present in the list of fields " << fields_
+ << exit(FatalIOError);
+ }
+
+
+ dict.lookup("trackForward") >> trackForward_;
+ dict.lookup("lifeTime") >> lifeTime_;
+ if (lifeTime_ < 1)
+ {
+ FatalErrorIn(":wallBoundedStreamLine::read(const dictionary&)")
+ << "Illegal value " << lifeTime_ << " for lifeTime"
+ << exit(FatalError);
+ }
+
+
+ interpolationScheme_ = dict.lookupOrDefault
+ (
+ "interpolationScheme",
+ interpolationCellPoint<scalar>::typeName
+ );
+
+ //Info<< typeName << " using interpolation " << interpolationScheme_
+ // << endl;
+
+ cloudName_ = dict.lookupOrDefault<word>
+ (
+ "cloudName",
+ "wallBoundedStreamLine"
+ );
+ dict.lookup("seedSampleSet") >> seedSet_;
+
+ const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+ const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
+
+ sampledSetPtr_ = sampledSet::New
+ (
+ seedSet_,
+ mesh,
+ meshSearchMeshObject::New(mesh),
+ coeffsDict
+ );
+ coeffsDict.lookup("axis") >> sampledSetAxis_;
+
+ scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
+ vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
+
+
+ // Make sure that the mesh is trackable
+ if (debug)
+ {
+ // 1. positive volume decomposition tets
+ faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
+ if
+ (
+ polyMeshTetDecomposition::checkFaceTets
+ (
+ mesh,
+ polyMeshTetDecomposition::minTetQuality,
+ true,
+ &faces
+ )
+ )
+ {
+ label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+ FatalErrorIn("wallBoundedStreamLine::track()")
+ << "Found " << nFaces
+ <<" faces with low quality or negative volume "
+ << "decomposition tets. Writing to faceSet " << faces.name()
+ << exit(FatalError);
+ }
+
+ // 2. all edges on a cell having two faces
+ EdgeMap<label> numFacesPerEdge;
+ forAll(mesh.cells(), cellI)
+ {
+ const cell& cFaces = mesh.cells()[cellI];
+
+ numFacesPerEdge.clear();
+
+ forAll(cFaces, cFaceI)
+ {
+ label faceI = cFaces[cFaceI];
+ const face& f = mesh.faces()[faceI];
+ forAll(f, fp)
+ {
+ const edge e(f[fp], f.nextLabel(fp));
+ EdgeMap<label>::iterator eFnd =
+ numFacesPerEdge.find(e);
+ if (eFnd != numFacesPerEdge.end())
+ {
+ eFnd()++;
+ }
+ else
+ {
+ numFacesPerEdge.insert(e, 1);
+ }
+ }
+ }
+
+ forAllConstIter(EdgeMap<label>, numFacesPerEdge, iter)
+ {
+ if (iter() != 2)
+ {
+ FatalErrorIn
+ (
+ "wallBoundedStreamLine::read(..)"
+ ) << "problem cell:" << cellI
+ << abort(FatalError);
+ }
+ }
+ }
+ }
+ }
+}
+
+
+void Foam::wallBoundedStreamLine::execute()
+{}
+
+
+void Foam::wallBoundedStreamLine::end()
+{}
+
+
+void Foam::wallBoundedStreamLine::write()
+{
+ if (active_)
+ {
+ const Time& runTime = obr_.time();
+ const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+
+ // Do all injection and tracking
+ track();
+
+
+ if (Pstream::parRun())
+ {
+ // Append slave tracks to master ones
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ globalIndex globalTrackIDs(allTracks_.size());
+
+ // Construct a distribution map to pull all to the master.
+ labelListList sendMap(Pstream::nProcs());
+ labelListList recvMap(Pstream::nProcs());
+
+ if (Pstream::master())
+ {
+ // Master: receive all. My own first, then consecutive
+ // processors.
+ label trackI = 0;
+
+ forAll(recvMap, procI)
+ {
+ labelList& fromProc = recvMap[procI];
+ fromProc.setSize(globalTrackIDs.localSize(procI));
+ forAll(fromProc, i)
+ {
+ fromProc[i] = trackI++;
+ }
+ }
+ }
+
+ labelList& toMaster = sendMap[0];
+ toMaster.setSize(globalTrackIDs.localSize());
+ forAll(toMaster, i)
+ {
+ toMaster[i] = i;
+ }
+
+ const mapDistribute distMap
+ (
+ globalTrackIDs.size(),
+ sendMap.xfer(),
+ recvMap.xfer()
+ );
+
+
+ // Distribute the track positions. Note: use scheduled comms
+ // to prevent buffering.
+ mapDistribute::distribute
+ (
+ Pstream::scheduled,
+ distMap.schedule(),
+ distMap.constructSize(),
+ distMap.subMap(),
+ distMap.constructMap(),
+ allTracks_
+ );
+
+ // Distribute the scalars
+ forAll(allScalars_, scalarI)
+ {
+ mapDistribute::distribute
+ (
+ Pstream::scheduled,
+ distMap.schedule(),
+ distMap.constructSize(),
+ distMap.subMap(),
+ distMap.constructMap(),
+ allScalars_[scalarI]
+ );
+ }
+ // Distribute the vectors
+ forAll(allVectors_, vectorI)
+ {
+ mapDistribute::distribute
+ (
+ Pstream::scheduled,
+ distMap.schedule(),
+ distMap.constructSize(),
+ distMap.subMap(),
+ distMap.constructMap(),
+ allVectors_[vectorI]
+ );
+ }
+ }
+
+
+ label n = 0;
+ forAll(allTracks_, trackI)
+ {
+ n += allTracks_[trackI].size();
+ }
+
+ Info<< "Tracks:" << allTracks_.size()
+ << " total samples:" << n << endl;
+
+
+ // Massage into form suitable for writers
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ if (Pstream::master() && allTracks_.size())
+ {
+ // Make output directory
+
+ fileName vtkPath
+ (
+ Pstream::parRun()
+ ? runTime.path()/".."/"sets"/name()
+ : runTime.path()/"sets"/name()
+ );
+ if (mesh.name() != fvMesh::defaultRegion)
+ {
+ vtkPath = vtkPath/mesh.name();
+ }
+ vtkPath = vtkPath/mesh.time().timeName();
+
+ mkDir(vtkPath);
+
+ // Convert track positions
+
+ PtrList<coordSet> tracks(allTracks_.size());
+ forAll(allTracks_, trackI)
+ {
+ tracks.set
+ (
+ trackI,
+ new coordSet
+ (
+ "track" + Foam::name(trackI),
+ sampledSetAxis_ //"xyz"
+ )
+ );
+ tracks[trackI].transfer(allTracks_[trackI]);
+ }
+
+ // Convert scalar values
+
+ if (allScalars_.size() > 0)
+ {
+ List<List<scalarField> > scalarValues(allScalars_.size());
+
+ forAll(allScalars_, scalarI)
+ {
+ DynamicList<scalarList>& allTrackVals =
+ allScalars_[scalarI];
+ scalarValues[scalarI].setSize(allTrackVals.size());
+
+ forAll(allTrackVals, trackI)
+ {
+ scalarList& trackVals = allTrackVals[trackI];
+ scalarValues[scalarI][trackI].transfer(trackVals);
+ }
+ }
+
+ fileName vtkFile
+ (
+ vtkPath
+ / scalarFormatterPtr_().getFileName
+ (
+ tracks[0],
+ scalarNames_
+ )
+ );
+
+ Info<< "Writing data to " << vtkFile.path() << endl;
+
+ scalarFormatterPtr_().write
+ (
+ true, // writeTracks
+ tracks,
+ scalarNames_,
+ scalarValues,
+ OFstream(vtkFile)()
+ );
+ }
+
+ // Convert vector values
+
+ if (allVectors_.size() > 0)
+ {
+ List<List<vectorField> > vectorValues(allVectors_.size());
+
+ forAll(allVectors_, vectorI)
+ {
+ DynamicList<vectorList>& allTrackVals =
+ allVectors_[vectorI];
+ vectorValues[vectorI].setSize(allTrackVals.size());
+
+ forAll(allTrackVals, trackI)
+ {
+ vectorList& trackVals = allTrackVals[trackI];
+ vectorValues[vectorI][trackI].transfer(trackVals);
+ }
+ }
+
+ fileName vtkFile
+ (
+ vtkPath
+ / vectorFormatterPtr_().getFileName
+ (
+ tracks[0],
+ vectorNames_
+ )
+ );
+
+ //Info<< "Writing vector data to " << vtkFile << endl;
+
+ vectorFormatterPtr_().write
+ (
+ true, // writeTracks
+ tracks,
+ vectorNames_,
+ vectorValues,
+ OFstream(vtkFile)()
+ );
+ }
+ }
+ }
+}
+
+
+void Foam::wallBoundedStreamLine::updateMesh(const mapPolyMesh&)
+{
+ read(dict_);
+}
+
+
+void Foam::wallBoundedStreamLine::movePoints(const pointField&)
+{
+ // Moving mesh affects the search tree
+ read(dict_);
+}
+
+
+//void Foam::wallBoundedStreamLine::readUpdate
+//(const polyMesh::readUpdateState state)
+//{
+// if (state != UNCHANGED)
+// {
+// read(dict_);
+// }
+//}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
new file mode 100644
index 0000000000000000000000000000000000000000..6a39d25b37da1948004ddf1973c8cc1cb78b2b3e
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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::wallBoundedStreamLine
+
+Description
+ Generation of streamlines. Samples along track of passive particle.
+
+SourceFiles
+ wallBoundedStreamLine.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLine_H
+#define wallBoundedStreamLine_H
+
+#include "volFieldsFwd.H"
+#include "pointFieldFwd.H"
+#include "Switch.H"
+#include "DynamicList.H"
+#include "scalarList.H"
+#include "vectorList.H"
+#include "polyMesh.H"
+#include "writer.H"
+#include "indirectPrimitivePatch.H"
+#include "tetIndices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+class meshSearch;
+class sampledSet;
+
+/*---------------------------------------------------------------------------*\
+ Class wallBoundedStreamLine Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLine
+{
+ // Private data
+
+ //- Input dictionary
+ dictionary dict_;
+
+ //- Name of this set of field averages.
+ word name_;
+
+ //- Database this class is registered to
+ const objectRegistry& obr_;
+
+ //- Load fields from files (not from objectRegistry)
+ bool loadFromFiles_;
+
+ //- On/off switch
+ bool active_;
+
+
+ //- List of fields to sample
+ wordList fields_;
+
+ //- Field to transport particle with
+ word UName_;
+
+ //- Interpolation scheme to use
+ word interpolationScheme_;
+
+ //- Whether to use +u or -u
+ bool trackForward_;
+
+ //- Maximum lifetime (= number of cells) of particle
+ label lifeTime_;
+
+ //- Optional specified name of particles
+ word cloudName_;
+
+ //- Type of seed
+ word seedSet_;
+
+ //- Names of scalar fields
+ wordList scalarNames_;
+
+ //- Names of vector fields
+ wordList vectorNames_;
+
+
+ // Demand driven
+
+ //- Mesh searching enigne
+ autoPtr<meshSearch> meshSearchPtr_;
+
+ //- Seed set engine
+ autoPtr<sampledSet> sampledSetPtr_;
+
+ //- Axis of the sampled points to output
+ word sampledSetAxis_;
+
+ //- File output writer
+ autoPtr<writer<scalar> > scalarFormatterPtr_;
+
+ autoPtr<writer<vector> > vectorFormatterPtr_;
+
+
+ // Generated data
+
+ //- All tracks. Per particle the points it passed through
+ DynamicList<List<point> > allTracks_;
+
+ //- Per scalarField, per particle, the sampled value.
+ List<DynamicList<scalarList> > allScalars_;
+
+ //- Per scalarField, per particle, the sampled value.
+ List<DynamicList<vectorList> > allVectors_;
+
+
+ //- Construct patch out of all wall patch faces
+ autoPtr<indirectPrimitivePatch> wallPatch() const;
+
+ //- Find wall tet on cell
+ tetIndices findNearestTet
+ (
+ const PackedBoolList& isWallPatch,
+ const point& seedPt,
+ const label cellI
+ ) const;
+
+ //- Do all seeding and tracking
+ void track();
+
+ //- Disallow default bitwise copy construct
+ wallBoundedStreamLine(const wallBoundedStreamLine&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const wallBoundedStreamLine&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("wallBoundedStreamLine");
+
+
+ // Constructors
+
+ //- Construct for given objectRegistry and dictionary.
+ // Allow the possibility to load fields from files
+ wallBoundedStreamLine
+ (
+ const word& name,
+ const objectRegistry&,
+ const dictionary&,
+ const bool loadFromFiles = false
+ );
+
+
+ //- Destructor
+ virtual ~wallBoundedStreamLine();
+
+
+ // Member Functions
+
+ //- Return name of the set of field averages
+ virtual const word& name() const
+ {
+ return name_;
+ }
+
+ //- Read the field average data
+ virtual void read(const dictionary&);
+
+ //- Execute the averaging
+ virtual void execute();
+
+ //- Execute the averaging at the final time-loop, currently does nothing
+ virtual void end();
+
+ //- Calculate the field average data and write
+ virtual void write();
+
+ //- Update for changes of mesh
+ virtual void updateMesh(const mapPolyMesh&);
+
+ //- Update for mesh point-motion
+ virtual void movePoints(const pointField&);
+
+ ////- Update for changes of mesh due to readUpdate
+ //virtual void readUpdate(const polyMesh::readUpdateState state);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
new file mode 100644
index 0000000000000000000000000000000000000000..eaa1cc7eaebbe0a7951ce5a1918021b7a51c21c2
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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 "wallBoundedStreamLineFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineNamedTemplateTypeNameAndDebug(wallBoundedStreamLineFunctionObject, 0);
+
+ addToRunTimeSelectionTable
+ (
+ functionObject,
+ wallBoundedStreamLineFunctionObject,
+ dictionary
+ );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H
new file mode 100644
index 0000000000000000000000000000000000000000..0ce95cc592a86b962edffd38271c508ab6795ed0
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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::wallBoundedStreamLineFunctionObject
+
+Description
+ FunctionObject wrapper around wallBoundedStreamLines
+ to allow them to be created via
+ the functions entry within controlDict.
+
+SourceFiles
+ wallBoundedStreamLineFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLineFunctionObject_H
+#define wallBoundedStreamLineFunctionObject_H
+
+#include "wallBoundedStreamLine.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ typedef OutputFilterFunctionObject<wallBoundedStreamLine>
+ wallBoundedStreamLineFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
new file mode 100644
index 0000000000000000000000000000000000000000..cd7605a56edd2c0b27eb172a48a9e0e060ec70cc
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -0,0 +1,1324 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*----------------------------------------------------------------------------*/
+
+#include "wallBoundedStreamLineParticle.H"
+#include "vectorFieldIOField.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+// defineParticleTypeNameAndDebug(wallBoundedStreamLineParticle, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+//// Check position is inside tet
+//void Foam::wallBoundedStreamLineParticle::checkInside() const
+//{
+// const tetIndices ti(currentTetIndices());
+// const tetPointRef tpr(ti.tet(mesh_));
+// if (!tpr.inside(position()))
+// {
+// FatalErrorIn("wallBoundedStreamLineParticle::checkInside(..)")
+// << "Particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "is not inside " << tpr
+// << abort(FatalError);
+// }
+//}
+//
+//
+//void Foam::wallBoundedStreamLineParticle::checkOnEdge() const
+//{
+// // Check that edge (as indicated by meshEdgeStart_, diagEdge_) is
+// // indeed one that contains the position.
+// const edge e = currentEdge();
+//
+// linePointRef ln(e.line(mesh_.points()));
+//
+// pointHit ph(ln.nearestDist(position()));
+//
+// if (ph.distance() > 1E-6)
+// {
+// FatalErrorIn
+// (
+// "wallBoundedStreamLineParticle::checkOnEdge()"
+// ) << "Problem :"
+// << " particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "edge:" << e
+// << " at:" << ln
+// << " distance:" << ph.distance()
+// << abort(FatalError);
+// }
+//}
+//
+//
+//void Foam::wallBoundedStreamLineParticle::checkOnTriangle(const point& p)
+//const
+//{
+// const triFace tri(currentTetIndices().faceTriIs(mesh_));
+// pointHit ph = tri.nearestPoint(p, mesh_.points());
+// if (ph.distance() > 1e-9)
+// {
+// FatalErrorIn
+// (
+// "wallBoundedStreamLineParticle::checkOnTriangle(const point&)"
+// ) << "Problem :"
+// << " particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "point:" << p
+// << " distance:" << ph.distance()
+// << abort(FatalError);
+// }
+//}
+
+
+// Construct the edge the particle is on (according to meshEdgeStart_,
+// diagEdge_)
+Foam::edge Foam::wallBoundedStreamLineParticle::currentEdge() const
+{
+ if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
+ {
+ FatalErrorIn("wallBoundedStreamLineParticle::currentEdge() const")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "cannot both be on a mesh edge and a face-diagonal edge."
+ << " meshEdgeStart_:" << meshEdgeStart_
+ << " diagEdge_:" << diagEdge_
+ << abort(FatalError);
+ }
+
+ const Foam::face& f = mesh_.faces()[tetFace()];
+
+ if (meshEdgeStart_ != -1)
+ {
+ return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ }
+ else
+ {
+ label faceBasePtI = mesh_.tetBasePtIs()[tetFace()];
+ label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+ return edge(f[faceBasePtI], f[diagPtI]);
+ }
+}
+
+
+void Foam::wallBoundedStreamLineParticle::crossEdgeConnectedFace
+(
+ const edge& meshEdge
+)
+{
+ //label oldFaceI = tetFace();
+
+ // Update tetFace, tetPt
+ particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
+
+ // Update face to be same as tracking one
+ face() = tetFace();
+
+
+ // And adapt meshEdgeStart_.
+ const Foam::face& f = mesh_.faces()[tetFace()];
+ label fp = findIndex(f, meshEdge[0]);
+
+ if (f.nextLabel(fp) == meshEdge[1])
+ {
+ meshEdgeStart_ = fp;
+ }
+ else
+ {
+ label fpMin1 = f.rcIndex(fp);
+
+ if (f[fpMin1] == meshEdge[1])
+ {
+ meshEdgeStart_ = fpMin1;
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "wallBoundedStreamLineParticle::crossEdgeConnectedFace"
+ "(const edge&)"
+ ) << "Problem :"
+ << " particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "face:" << tetFace()
+ << " verts:" << f
+ << " meshEdge:" << meshEdge
+ << abort(FatalError);
+ }
+ }
+
+ diagEdge_ = -1;
+
+ //Pout<< " crossed meshEdge "
+ // << meshEdge.line(mesh().points())
+ // << " from face:" << oldFaceI
+ // << " to face:" << tetFace() << endl;
+
+
+ // Check that still on same mesh edge
+
+ const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ if (eNew != meshEdge)
+ {
+ FatalErrorIn
+ (
+ "wallBoundedStreamLineParticle::crossEdgeConnectedFace"
+ "(const edge&)"
+ ) << "Problem" << abort(FatalError);
+ }
+
+
+ // Check that edge (as indicated by meshEdgeStart_) is indeed one that
+ // contains the position.
+ //checkOnEdge();
+}
+
+
+void Foam::wallBoundedStreamLineParticle::crossDiagonalEdge()
+{
+ if (diagEdge_ == -1)
+ {
+ FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "not on a diagonal edge" << abort(FatalError);
+ }
+ if (meshEdgeStart_ != -1)
+ {
+ FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
+ }
+
+ //label oldTetPt = tetPt();
+
+ const Foam::face& f = mesh_.faces()[tetFace()];
+
+ // tetPtI starts from 1, goes up to f.size()-2
+
+ if (tetPt() == diagEdge_)
+ {
+ tetPt() = f.rcIndex(tetPt());
+ }
+ else
+ {
+ label nextTetPt = f.fcIndex(tetPt());
+ if (diagEdge_ == nextTetPt)
+ {
+ tetPt() = nextTetPt;
+ }
+ else
+ {
+ FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "tetPt:" << tetPt()
+ << " diagEdge:" << diagEdge_ << abort(FatalError);
+ }
+ }
+
+ meshEdgeStart_ = -1;
+
+ //Pout<< " crossed diagonalEdge "
+ // << currentEdge().line(mesh().points())
+ // << " from tetPt:" << oldTetPt
+ // << " to tetPt:" << tetPt() << endl;
+}
+
+
+//- Track through a single triangle.
+// Gets passed tet+triangle the particle is in. Updates position() but nothing
+// else. Returns the triangle edge the particle is now on.
+Foam::scalar Foam::wallBoundedStreamLineParticle::trackFaceTri
+(
+ const vector& endPosition,
+ label& minEdgeI
+)
+{
+ // Track p from position to endPosition
+ const triFace tri(currentTetIndices().faceTriIs(mesh_));
+ vector n = tri.normal(mesh_.points());
+ //if (mag(n) < sqr(SMALL))
+ //{
+ // FatalErrorIn("wallBoundedStreamLineParticle::trackFaceTri(..)")
+ // << "Small triangle." //<< static_cast<const particle&>(*this)
+ // << info()
+ // << "n:" << n
+ // << abort(FatalError);
+ //}
+ n /= mag(n)+VSMALL;
+
+ // Check which edge intersects the trajectory.
+ // Project trajectory onto triangle
+ minEdgeI = -1;
+ scalar minS = 1; // end position
+
+ //const point oldPosition(position());
+
+
+ edge currentE(-1, -1);
+ if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ {
+ currentE = currentEdge();
+ }
+
+ // Determine path along line position+s*d to see where intersections
+ // are.
+
+ forAll(tri, i)
+ {
+ label j = tri.fcIndex(i);
+
+ const point& pt0 = mesh_.points()[tri[i]];
+ const point& pt1 = mesh_.points()[tri[j]];
+
+ if (edge(tri[i], tri[j]) == currentE)
+ {
+ // Do not check particle is on
+ continue;
+ }
+
+ // Outwards pointing normal
+ vector edgeNormal = (pt1-pt0)^n;
+
+ //if (mag(edgeNormal) < SMALL)
+ //{
+ // FatalErrorIn("wallBoundedStreamLineParticle::trackFaceTri(..)")
+ // << "Edge not perpendicular to triangle."
+ // //<< static_cast<const particle&>(*this)
+ // << info()
+ // << "triangle n:" << n
+ // << " edgeNormal:" << edgeNormal
+ // << " on tri:" << tri
+ // << " at:" << pt0
+ // << " at:" << pt1
+ // << abort(FatalError);
+ //}
+
+
+ edgeNormal /= mag(edgeNormal)+VSMALL;
+
+ // Determine whether position and end point on either side of edge.
+ scalar sEnd = (endPosition-pt0)&edgeNormal;
+ if (sEnd >= 0)
+ {
+ // endPos is outside triangle. position() should always be
+ // inside.
+ scalar sStart = (position()-pt0)&edgeNormal;
+ if (mag(sEnd-sStart) > VSMALL)
+ {
+ scalar s = sStart/(sStart-sEnd);
+
+ if (s >= 0 && s < minS)
+ {
+ minS = s;
+ minEdgeI = i;
+ }
+ }
+ }
+ }
+
+ if (minEdgeI != -1)
+ {
+ position() += minS*(endPosition-position());
+ }
+ else
+ {
+ // Did not hit any edge so tracked to the end position. Set position
+ // without any calculation to avoid truncation errors.
+ position() = endPosition;
+ minS = 1.0;
+ }
+
+ // Project position() back onto plane of triangle
+ const point& triPt = mesh_.points()[tri[0]];
+ position() -= ((position()-triPt)&n)*n;
+
+
+ //Pout<< " tracked from:" << oldPosition << " to:" << position()
+ // << " projectedEnd:" << endPosition
+ // << " at s:" << minS << endl;
+ //if (minEdgeI != -1)
+ //{
+ // Pout<< " on edge:" << minEdgeI
+ // << " on edge:"
+ // << mesh_.points()[tri[minEdgeI]]
+ // << mesh_.points()[tri[tri.fcIndex(minEdgeI)]]
+ // << endl;
+ //}
+
+ return minS;
+}
+
+
+// See if the current triangle has got a point on the
+// correct side of the edge.
+bool Foam::wallBoundedStreamLineParticle::isTriAlongTrack
+(
+ const point& endPosition
+) const
+{
+ const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
+ const edge currentE = currentEdge();
+
+ //if (debug)
+ {
+ if
+ (
+ currentE[0] == currentE[1]
+ || findIndex(triVerts, currentE[0]) == -1
+ || findIndex(triVerts, currentE[1]) == -1
+ )
+ {
+ FatalErrorIn
+ (
+ "wallBoundedStreamLineParticle::isTriAlongTrack"
+ "(const point&)"
+ ) << "Edge " << currentE << " not on triangle " << triVerts
+ << info()
+ << abort(FatalError);
+ }
+ }
+
+
+ const vector dir = endPosition-position();
+
+ // Get normal of currentE
+ vector n = triVerts.normal(mesh_.points());
+ n /= mag(n);
+
+ forAll(triVerts, i)
+ {
+ label j = triVerts.fcIndex(i);
+ const point& pt0 = mesh_.points()[triVerts[i]];
+ const point& pt1 = mesh_.points()[triVerts[j]];
+
+ if (edge(triVerts[i], triVerts[j]) == currentE)
+ {
+ vector edgeNormal = (pt1-pt0)^n;
+ return (dir&edgeNormal) < 0;
+ }
+ }
+
+ FatalErrorIn
+ (
+ "wallBoundedStreamLineParticle::isTriAlongTrack"
+ "(const point&)"
+ ) << "Problem" << abort(FatalError);
+
+ return false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::patchInteraction
+(
+ trackingData& td,
+ const scalar trackFraction
+)
+{
+ typedef trackingData::cloudType cloudType;
+ typedef cloudType::particleType particleType;
+
+ particleType& p = static_cast<particleType&>(*this);
+ p.hitFace(td);
+
+ if (!internalFace(faceI_))
+ {
+ label origFaceI = faceI_;
+ label patchI = patch(faceI_);
+
+ // No action taken for tetPtI_ for tetFaceI_ here, handled by
+ // patch interaction call or later during processor transfer.
+
+
+ // Dummy tet indices. What to do here?
+ tetIndices faceHitTetIs;
+
+ if
+ (
+ !p.hitPatch
+ (
+ mesh_.boundaryMesh()[patchI],
+ td,
+ patchI,
+ trackFraction,
+ faceHitTetIs
+ )
+ )
+ {
+ // Did patch interaction model switch patches?
+ // Note: recalculate meshEdgeStart_, diagEdge_!
+ if (faceI_ != origFaceI)
+ {
+ patchI = patch(faceI_);
+ }
+
+ const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+
+ if (isA<wedgePolyPatch>(patch))
+ {
+ p.hitWedgePatch
+ (
+ static_cast<const wedgePolyPatch&>(patch), td
+ );
+ }
+ else if (isA<symmetryPolyPatch>(patch))
+ {
+ p.hitSymmetryPatch
+ (
+ static_cast<const symmetryPolyPatch&>(patch), td
+ );
+ }
+ else if (isA<cyclicPolyPatch>(patch))
+ {
+ p.hitCyclicPatch
+ (
+ static_cast<const cyclicPolyPatch&>(patch), td
+ );
+ }
+ else if (isA<processorPolyPatch>(patch))
+ {
+ p.hitProcessorPatch
+ (
+ static_cast<const processorPolyPatch&>(patch), td
+ );
+ }
+ else if (isA<wallPolyPatch>(patch))
+ {
+ p.hitWallPatch
+ (
+ static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
+ );
+ }
+ else
+ {
+ p.hitPatch(patch, td);
+ }
+ }
+ }
+}
+
+//- 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.
+Foam::scalar Foam::wallBoundedStreamLineParticle::trackToEdge
+(
+ trackingData& td,
+ const vector& endPosition
+)
+{
+ // Are we on a track face? If not we do a topological walk.
+
+ // Particle:
+ // - cell_ always set
+ // - tetFace_, tetPt_ always set (these identify tet particle is in)
+ // - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
+
+ //checkInside();
+ //checkOnTriangle(position());
+ //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ //{
+ // checkOnEdge();
+ //}
+
+ scalar trackFraction = 0.0;
+
+ if (!td.isWallPatch_[tetFace()])
+ {
+ // Don't track across face. Just walk in cell. Particle is on
+ // mesh edge (as indicated by meshEdgeStart_).
+ const edge meshEdge(currentEdge());
+
+ // If internal face check whether to go to neighbour cell or just
+ // check to the other internal tet on the edge.
+ if (mesh_.isInternalFace(tetFace()))
+ {
+ label nbrCellI =
+ (
+ cellI_ == mesh_.faceOwner()[faceI_]
+ ? mesh_.faceNeighbour()[faceI_]
+ : mesh_.faceOwner()[faceI_]
+ );
+ // Check angle to nbrCell tet. Is it in the direction of the
+ // endposition? I.e. since volume of nbr tet is positive the
+ // tracking direction should be into the tet.
+ tetIndices nbrTi(nbrCellI, tetFaceI_, tetPtI_, mesh_);
+ if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
+ {
+ // Change into nbrCell. No need to change tetFace, tetPt.
+ //Pout<< " crossed from cell:" << cellI_
+ // << " into " << nbrCellI << endl;
+ cellI_ = nbrCellI;
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Walk to other face on edge. Changes tetFace, tetPt but not
+ // cell.
+ crossEdgeConnectedFace(meshEdge);
+ patchInteraction(td, trackFraction);
+ }
+ }
+ else
+ {
+ // Walk to other face on edge. This might give loop since
+ // particle should have been removed?
+ crossEdgeConnectedFace(meshEdge);
+ patchInteraction(td, trackFraction);
+ }
+ }
+ else
+ {
+ // We're inside a tet on the wall. Check if the current tet is
+ // the one to cross. If not we cross into the neighbouring triangle.
+
+ if (mesh_.isInternalFace(tetFace()))
+ {
+ FatalErrorIn
+ (
+ "wallBoundedStreamLineParticle::trackToEdge"
+ "(trackingData&, const vector&)"
+ ) << "Can only track on boundary faces."
+ << " Face:" << tetFace()
+ << " at:" << mesh_.faceCentres()[tetFace()]
+ << abort(FatalError);
+ }
+
+ point projectedEndPosition = endPosition;
+ // Remove normal component
+ {
+ const triFace tri(currentTetIndices().faceTriIs(mesh_));
+ vector n = tri.normal(mesh_.points());
+ n /= mag(n);
+ const point& basePt = mesh_.points()[tri[0]];
+ projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
+ }
+
+
+ bool doTrack = false;
+ if (meshEdgeStart_ == -1 && diagEdge_ == -1)
+ {
+ // We're starting and not yet on an edge.
+ doTrack = true;
+ }
+ else
+ {
+ // See if the current triangle has got a point on the
+ // correct side of the edge.
+ doTrack = isTriAlongTrack(projectedEndPosition);
+ }
+
+
+ if (doTrack)
+ {
+ // Track across triangle. Return triangle edge crossed.
+ label triEdgeI = -1;
+ trackFraction = trackFaceTri(projectedEndPosition, triEdgeI);
+
+ if (triEdgeI == -1)
+ {
+ // Reached endpoint
+ //checkInside();
+ return trackFraction;
+ }
+
+ const tetIndices ti(currentTetIndices());
+
+ // Triangle (faceTriIs) gets constructed from
+ // f[faceBasePtI_],
+ // f[facePtAI_],
+ // f[facePtBI_]
+ //
+ // So edge indices are:
+ // 0 : edge between faceBasePtI_ and facePtAI_
+ // 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
+ // 2 : edge between facePtBI_ and faceBasePtI_
+
+ const Foam::face& f = mesh_.faces()[ti.face()];
+ const label fp0 = ti.faceBasePt();
+
+ if (triEdgeI == 0)
+ {
+ if (ti.facePtA() == f.fcIndex(fp0))
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = fp0;
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else if (ti.facePtA() == f.rcIndex(fp0))
+ {
+ //Note: should not happen since boundary face so owner
+ //Pout<< "Real edge." << endl;
+ FatalErrorIn("shold not happend") << info()
+ << abort(FatalError);
+
+ diagEdge_ = -1;
+ meshEdgeStart_ = f.rcIndex(fp0);
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Get index of triangle on other side of edge.
+ diagEdge_ = ti.facePtA()-fp0;
+ if (diagEdge_ < 0)
+ {
+ diagEdge_ += f.size();
+ }
+ meshEdgeStart_ = -1;
+ //checkOnEdge();
+ crossDiagonalEdge();
+ }
+ }
+ else if (triEdgeI == 1)
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = ti.facePtA();
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else // if (triEdgeI == 2)
+ {
+ if (ti.facePtB() == f.rcIndex(fp0))
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = ti.facePtB();
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else if (ti.facePtB() == f.fcIndex(fp0))
+ {
+ //Note: should not happen since boundary face so owner
+ //Pout<< "Real edge." << endl;
+ FatalErrorIn("shold not happend") << info()
+ << abort(FatalError);
+
+ diagEdge_ = -1;
+ meshEdgeStart_ = fp0;
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ //Pout<< "Triangle edge." << endl;
+ // Get index of triangle on other side of edge.
+ diagEdge_ = ti.facePtB()-fp0;
+ if (diagEdge_ < 0)
+ {
+ diagEdge_ += f.size();
+ }
+ meshEdgeStart_ = -1;
+ //checkOnEdge();
+ crossDiagonalEdge();
+ }
+ }
+ }
+ else
+ {
+ // Current tet is not the right one. Check the neighbour tet.
+
+ if (meshEdgeStart_ != -1)
+ {
+ // Particle is on mesh edge so change into other face on cell
+ crossEdgeConnectedFace(currentEdge());
+ //checkOnEdge();
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Particle is on diagonal edge so change into the other
+ // triangle.
+ crossDiagonalEdge();
+ //checkOnEdge();
+ }
+ }
+ }
+
+ //checkInside();
+
+ return trackFraction;
+}
+
+
+Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
+(
+ const trackingData& td,
+ const point& position,
+ const label cellI,
+ const label faceI
+)
+{
+ if (cellI == -1)
+ {
+ FatalErrorIn("wallBoundedStreamLineParticle::interpolateFields(..)")
+ << "Cell:" << cellI << abort(FatalError);
+ }
+
+ const tetIndices ti = currentTetIndices();
+
+
+ sampledScalars_.setSize(td.vsInterp_.size());
+ forAll(td.vsInterp_, scalarI)
+ {
+ sampledScalars_[scalarI].append
+ (
+ td.vsInterp_[scalarI].interpolate
+ (
+ position,
+ ti, //cellI,
+ faceI
+ )
+ );
+ }
+
+ sampledVectors_.setSize(td.vvInterp_.size());
+ forAll(td.vvInterp_, vectorI)
+ {
+ sampledVectors_[vectorI].append
+ (
+ td.vvInterp_[vectorI].interpolate
+ (
+ position,
+ ti, //cellI,
+ faceI
+ )
+ );
+ }
+
+ const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
+
+ return U.last();
+}
+
+
+Foam::vector Foam::wallBoundedStreamLineParticle::sample
+(
+ trackingData& td
+)
+{
+ sampledPositions_.append(position());
+ vector U = interpolateFields(td, position(), cell(), tetFace());
+
+ if (!td.trackForward_)
+ {
+ U = -U;
+ }
+
+ scalar magU = mag(U);
+
+ if (magU < SMALL)
+ {
+ // Stagnant particle. Might as well stop
+ lifeTime_ = 0;
+ }
+
+ U /= magU;
+
+ return U;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+ const polyMesh& mesh,
+ const vector& position,
+ const label cellI,
+ const label tetFaceI,
+ const label tetPtI,
+ const label meshEdgeStart,
+ const label diagEdge,
+ const label lifeTime
+)
+:
+ particle(mesh, position, cellI, tetFaceI, tetPtI),
+ meshEdgeStart_(meshEdgeStart),
+ diagEdge_(diagEdge),
+ lifeTime_(lifeTime)
+{
+ //checkInside();
+
+ //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ //{
+ // checkOnEdge();
+ //}
+
+ // Unfortunately have no access to trackdata so cannot check if particle
+ // is on a wallPatch or has an mesh edge set (either of which is
+ // a requirement).
+}
+
+
+Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+ const polyMesh& mesh,
+ Istream& is,
+ bool readFields
+)
+:
+ particle(mesh, is, readFields)
+{
+ if (readFields)
+ {
+ List<scalarList> sampledScalars;
+ List<vectorList> sampledVectors;
+
+ is >> meshEdgeStart_ >> diagEdge_ >> lifeTime_
+ >> sampledPositions_ >> sampledScalars >> sampledVectors;
+
+ sampledScalars_.setSize(sampledScalars.size());
+ forAll(sampledScalars, i)
+ {
+ sampledScalars_[i].transfer(sampledScalars[i]);
+ }
+ sampledVectors_.setSize(sampledVectors.size());
+ forAll(sampledVectors, i)
+ {
+ sampledVectors_[i].transfer(sampledVectors[i]);
+ }
+ }
+
+ // Check state of Istream
+ is.check
+ (
+ "wallBoundedStreamLineParticle::wallBoundedStreamLineParticle"
+ "(const Cloud<wallBoundedStreamLineParticle>&, Istream&, bool)"
+ );
+}
+
+
+Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+ const wallBoundedStreamLineParticle& p
+)
+:
+ particle(p),
+ meshEdgeStart_(p.meshEdgeStart_),
+ diagEdge_(p.diagEdge_),
+ lifeTime_(p.lifeTime_),
+ sampledPositions_(p.sampledPositions_),
+ sampledScalars_(p.sampledScalars_),
+ sampledVectors_(p.sampledVectors_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::wallBoundedStreamLineParticle::move
+(
+ trackingData& td,
+ const scalar trackTime
+)
+{
+ wallBoundedStreamLineParticle& p = static_cast
+ <
+ wallBoundedStreamLineParticle&
+ >(*this);
+
+
+ // Check position is inside tet
+ //checkInside();
+
+ td.switchProcessor = false;
+ td.keepParticle = true;
+
+ scalar tEnd = (1.0 - stepFraction())*trackTime;
+ scalar maxDt = mesh_.bounds().mag();
+
+ while
+ (
+ td.keepParticle
+ && !td.switchProcessor
+ && lifeTime_ > 0
+ )
+ {
+ // set the lagrangian time-step
+ scalar dt = maxDt;
+
+ --lifeTime_;
+
+ // Update sampled velocity and fields if position changed
+ vector U = sample(td);
+
+ // !user parameter!
+ if (dt < SMALL)
+ {
+ // Force removal
+ lifeTime_ = 0;
+ break;
+ }
+
+
+ scalar fraction = trackToEdge(td, position() + dt*U);
+ dt *= fraction;
+
+ tEnd -= dt;
+ stepFraction() = 1.0 - tEnd/trackTime;
+
+
+ if (tEnd <= ROOTVSMALL)
+ {
+ // Force removal
+ lifeTime_ = 0;
+ }
+
+ if
+ (
+ !td.keepParticle
+ || td.switchProcessor
+ || lifeTime_ == 0
+ )
+ {
+ break;
+ }
+ }
+
+
+ if (!td.keepParticle || lifeTime_ == 0)
+ {
+ if (lifeTime_ == 0)
+ {
+ if (debug)
+ {
+ Pout<< "wallBoundedStreamLineParticle :"
+ << " Removing stagnant particle:"
+ << p.position()
+ << " sampled positions:" << sampledPositions_.size()
+ << endl;
+ }
+ td.keepParticle = false;
+ }
+ else
+ {
+ // Normal exit. Store last position and fields
+ sample(td);
+
+ if (debug)
+ {
+ Pout<< "wallBoundedStreamLineParticle : Removing particle:"
+ << p.position()
+ << " sampled positions:" << sampledPositions_.size()
+ << endl;
+ }
+ }
+
+ // Transfer particle data into trackingData.
+ {
+ //td.allPositions_.append(sampledPositions_);
+ td.allPositions_.append(vectorList());
+ vectorList& top = td.allPositions_.last();
+ top.transfer(sampledPositions_);
+ }
+
+ forAll(sampledScalars_, i)
+ {
+ //td.allScalars_[i].append(sampledScalars_[i]);
+ td.allScalars_[i].append(scalarList());
+ scalarList& top = td.allScalars_[i].last();
+ top.transfer(sampledScalars_[i]);
+ }
+ forAll(sampledVectors_, i)
+ {
+ //td.allVectors_[i].append(sampledVectors_[i]);
+ td.allVectors_[i].append(vectorList());
+ vectorList& top = td.allVectors_[i].last();
+ top.transfer(sampledVectors_[i]);
+ }
+ }
+
+ return td.keepParticle;
+}
+
+
+bool Foam::wallBoundedStreamLineParticle::hitPatch
+(
+ const polyPatch&,
+ trackingData& td,
+ const label patchI,
+ const scalar trackFraction,
+ const tetIndices& tetIs
+)
+{
+ // Disable generic patch interaction
+ return false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitWedgePatch
+(
+ const wedgePolyPatch& pp,
+ trackingData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitSymmetryPatch
+(
+ const symmetryPolyPatch& pp,
+ trackingData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitCyclicPatch
+(
+ const cyclicPolyPatch& pp,
+ trackingData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitProcessorPatch
+(
+ const processorPolyPatch& pp,
+ trackingData& td
+)
+{
+ // Switch particle
+ td.switchProcessor = true;
+
+ // Adapt edgeStart_ for other side.
+ // E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
+ // on the other side between 2 and 3 so edgeStart_ should be
+ // f.size()-edgeStart_-1.
+
+ const Foam::face& f = mesh_.faces()[face()];
+
+ if (meshEdgeStart_ != -1)
+ {
+ meshEdgeStart_ = f.size()-meshEdgeStart_-1;
+ }
+ else
+ {
+ // diagEdge_ is relative to faceBasePt
+ diagEdge_ = f.size()-diagEdge_;
+ }
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitWallPatch
+(
+ const wallPolyPatch& wpp,
+ trackingData& td,
+ const tetIndices&
+)
+{}
+
+
+void Foam::wallBoundedStreamLineParticle::hitPatch
+(
+ const polyPatch& wpp,
+ trackingData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::readFields
+(
+ Cloud<wallBoundedStreamLineParticle>& c
+)
+{
+ if (!c.size())
+ {
+ return;
+ }
+
+ particle::readFields(c);
+
+ IOField<label> meshEdgeStart
+ (
+ c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
+ );
+
+ IOField<label> diagEdge
+ (
+ c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
+ );
+ c.checkFieldIOobject(c, diagEdge);
+
+ IOField<label> lifeTime
+ (
+ c.fieldIOobject("lifeTime", IOobject::MUST_READ)
+ );
+ c.checkFieldIOobject(c, lifeTime);
+
+ vectorFieldIOField sampledPositions
+ (
+ c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
+ );
+ c.checkFieldIOobject(c, sampledPositions);
+
+ label i = 0;
+ forAllIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
+ {
+ iter().meshEdgeStart_ = meshEdgeStart[i];
+ iter().diagEdge_ = diagEdge[i];
+ iter().lifeTime_ = lifeTime[i];
+ iter().sampledPositions_.transfer(sampledPositions[i]);
+ i++;
+ }
+}
+
+
+void Foam::wallBoundedStreamLineParticle::writeFields
+(
+ const Cloud<wallBoundedStreamLineParticle>& c
+)
+{
+ particle::writeFields(c);
+
+ label np = c.size();
+
+ IOField<label> meshEdgeStart
+ (
+ c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
+ np
+ );
+ IOField<label> diagEdge
+ (
+ c.fieldIOobject("diagEdge", IOobject::NO_READ),
+ np
+ );
+ IOField<label> lifeTime
+ (
+ c.fieldIOobject("lifeTime", IOobject::NO_READ),
+ np
+ );
+ vectorFieldIOField sampledPositions
+ (
+ c.fieldIOobject("sampledPositions", IOobject::NO_READ),
+ np
+ );
+
+ label i = 0;
+ forAllConstIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
+ {
+ meshEdgeStart[i] = iter().meshEdgeStart_;
+ diagEdge[i] = iter().diagEdge_;
+ lifeTime[i] = iter().lifeTime_;
+ sampledPositions[i] = iter().sampledPositions_;
+ i++;
+ }
+
+ meshEdgeStart.write();
+ diagEdge.write();
+ lifeTime.write();
+ sampledPositions.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const wallBoundedStreamLineParticle& p
+)
+{
+ os << static_cast<const particle&>(p)
+ << token::SPACE << p.meshEdgeStart_
+ << token::SPACE << p.diagEdge_
+ << token::SPACE << p.lifeTime_
+ << token::SPACE << p.sampledPositions_
+ << token::SPACE << p.sampledScalars_
+ << token::SPACE << p.sampledVectors_;
+
+ // Check state of Ostream
+ os.check
+ (
+ "Ostream& operator<<(Ostream&, const wallBoundedStreamLineParticle&)"
+ );
+
+ return os;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const InfoProxy<wallBoundedStreamLineParticle>& ip
+)
+{
+ const wallBoundedStreamLineParticle& p = ip.t_;
+
+ tetPointRef tpr(p.currentTet());
+
+ os << " " << static_cast<const particle&>(p) << nl
+ << " tet:" << nl;
+ os << " ";
+ meshTools::writeOBJ(os, tpr.a());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.b());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.c());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.d());
+ os << " l 1 2" << nl
+ << " l 1 3" << nl
+ << " l 1 4" << nl
+ << " l 2 3" << nl
+ << " l 2 4" << nl
+ << " l 3 4" << nl;
+ os << " ";
+ meshTools::writeOBJ(os, p.position());
+
+ return os;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
new file mode 100644
index 0000000000000000000000000000000000000000..c0e19dd090a55e559f54e9b9d1e0d2562a307e05
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
@@ -0,0 +1,368 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::wallBoundedStreamLineParticle
+
+Description
+ Particle class that samples fields as it passes through. Used in streamline
+ calculation.
+
+SourceFiles
+ wallBoundedStreamLineParticle.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLineParticle_H
+#define wallBoundedStreamLineParticle_H
+
+#include "particle.H"
+#include "autoPtr.H"
+#include "interpolation.H"
+#include "vectorList.H"
+#include "InfoProxy.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class wallBoundedStreamLineParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+ Class wallBoundedStreamLineParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLineParticle
+:
+ public particle
+{
+
+public:
+
+ //- Class used to pass tracking data to the trackToFace function
+ class trackingData
+ :
+ public particle::TrackingData<Cloud<wallBoundedStreamLineParticle> >
+ {
+
+ public:
+
+
+ const PtrList<interpolation<scalar> >& vsInterp_;
+ const PtrList<interpolation<vector> >& vvInterp_;
+ const label UIndex_;
+ const bool trackForward_;
+
+ const PackedBoolList& isWallPatch_;
+
+ DynamicList<vectorList>& allPositions_;
+ List<DynamicList<scalarList> >& allScalars_;
+ List<DynamicList<vectorList> >& allVectors_;
+
+
+ // Constructors
+
+ trackingData
+ (
+ Cloud<wallBoundedStreamLineParticle>& cloud,
+ const PtrList<interpolation<scalar> >& vsInterp,
+ const PtrList<interpolation<vector> >& vvInterp,
+ const label UIndex,
+ const bool trackForward,
+ const PackedBoolList& isWallPatch,
+
+ DynamicList<List<point> >& allPositions,
+ List<DynamicList<scalarList> >& allScalars,
+ List<DynamicList<vectorList> >& allVectors
+ )
+ :
+ particle::TrackingData<Cloud<wallBoundedStreamLineParticle> >
+ (
+ cloud
+ ),
+ vsInterp_(vsInterp),
+ vvInterp_(vvInterp),
+ UIndex_(UIndex),
+ trackForward_(trackForward),
+ isWallPatch_(isWallPatch),
+
+ allPositions_(allPositions),
+ allScalars_(allScalars),
+ allVectors_(allVectors)
+ {}
+ };
+
+
+private:
+
+ // Private data
+
+ //- Particle is on mesh edge:
+ // const face& f = mesh.faces()[tetFace()]
+ // const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ // Note that this real edge
+ // is also one of the edges of the face-triangle (from
+ // tetFace()+tetPt()).
+ label meshEdgeStart_;
+
+ //- Particle is on diagonal edge:
+ // const face& f = mesh.faces()[tetFace()]
+ // label faceBasePtI = mesh.tetBasePtIs()[faceI];
+ // label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+ // const edge e(f[faceBasePtI], f[diagPtI]);
+ label diagEdge_;
+
+ //- Lifetime of particle. Particle dies when reaches 0.
+ label lifeTime_;
+
+ //- sampled positions
+ DynamicList<point> sampledPositions_;
+
+ //- sampled scalars
+ List<DynamicList<scalar> > sampledScalars_;
+
+ //- sampled vectors
+ List<DynamicList<vector> > sampledVectors_;
+
+
+ // Private Member Functions
+
+ //- Construct current edge
+ edge currentEdge() const;
+
+ //- Check if inside current tet
+ //void checkInside() const;
+
+ //- Check if on current edge
+ //void checkOnEdge() const;
+
+ //- Check if point on triangle
+ //void checkOnTriangle(const point&) const;
+
+ //- Cross mesh edge into different face on same cell
+ void crossEdgeConnectedFace(const edge& meshEdge);
+
+ //- Cross diagonal edge into different triangle on same face,cell
+ void crossDiagonalEdge();
+
+ //- Track through single triangle
+ scalar trackFaceTri(const vector& endPosition, label& minEdgeI);
+
+ //- Do all patch interaction
+ void patchInteraction(trackingData& td, const scalar trackFraction);
+
+ //- Is current triangle in the track direction
+ bool isTriAlongTrack(const point& endPosition) const;
+
+ //- Equivalent of trackToFace
+ scalar trackToEdge
+ (
+ trackingData& td,
+ const vector& endPosition
+ );
+
+ //- Interpolate all quantities; return interpolated velocity.
+ vector interpolateFields
+ (
+ const trackingData&,
+ const point&,
+ const label cellI,
+ const label faceI
+ );
+
+ //- Interpolate all quantities, return updated velocity.
+ vector sample(trackingData& td);
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ wallBoundedStreamLineParticle
+ (
+ const polyMesh& c,
+ const vector& position,
+ const label cellI,
+ const label tetFaceI,
+ const label tetPtI,
+ const label meshEdgeStart,
+ const label diagEdge,
+ const label lifeTime
+ );
+
+ //- Construct from Istream
+ wallBoundedStreamLineParticle
+ (
+ const polyMesh& c,
+ Istream& is,
+ bool readFields = true
+ );
+
+ //- Construct copy
+ wallBoundedStreamLineParticle(const wallBoundedStreamLineParticle& p);
+
+ //- Construct and return a clone
+ autoPtr<particle> clone() const
+ {
+ return autoPtr<particle>(new wallBoundedStreamLineParticle(*this));
+ }
+
+ //- Factory class to read-construct particles used for
+ // parallel transfer
+ class iNew
+ {
+ const polyMesh& mesh_;
+
+ public:
+
+ iNew(const polyMesh& mesh)
+ :
+ mesh_(mesh)
+ {}
+
+ autoPtr<wallBoundedStreamLineParticle> operator()
+ (
+ Istream& is
+ ) const
+ {
+ return autoPtr<wallBoundedStreamLineParticle>
+ (
+ new wallBoundedStreamLineParticle(mesh_, is, true)
+ );
+ }
+ };
+
+
+ // Member Functions
+
+ // Tracking
+
+ //- Track all particles to their end point
+ bool move(trackingData&, const scalar trackTime);
+
+
+ //- Overridable function to handle the particle hitting a patch
+ // Executed before other patch-hitting functions
+ bool hitPatch
+ (
+ const polyPatch&,
+ trackingData& td,
+ const label patchI,
+ const scalar trackFraction,
+ const tetIndices& tetIs
+ );
+
+ //- Overridable function to handle the particle hitting a wedge
+ void hitWedgePatch
+ (
+ const wedgePolyPatch&,
+ trackingData& td
+ );
+
+ //- Overridable function to handle the particle hitting a
+ // symmetryPlane
+ void hitSymmetryPatch
+ (
+ const symmetryPolyPatch&,
+ trackingData& td
+ );
+
+ //- Overridable function to handle the particle hitting a cyclic
+ void hitCyclicPatch
+ (
+ const cyclicPolyPatch&,
+ trackingData& td
+ );
+
+ //- Overridable function to handle the particle hitting a
+ //- processorPatch
+ void hitProcessorPatch
+ (
+ const processorPolyPatch&,
+ trackingData& td
+ );
+
+ //- Overridable function to handle the particle hitting a wallPatch
+ void hitWallPatch
+ (
+ const wallPolyPatch&,
+ trackingData& td,
+ const tetIndices&
+ );
+
+ //- Overridable function to handle the particle hitting a polyPatch
+ void hitPatch
+ (
+ const polyPatch&,
+ trackingData& td
+ );
+
+
+ // Info
+
+ //- Return info proxy.
+ // Used to print particle information to a stream
+ InfoProxy<wallBoundedStreamLineParticle> info() const
+ {
+ return *this;
+ }
+
+
+ // I-O
+
+ //- Read
+ static void readFields(Cloud<wallBoundedStreamLineParticle>&);
+
+ //- Write
+ static void writeFields
+ (
+ const Cloud<wallBoundedStreamLineParticle>&
+ );
+
+
+ // Ostream Operator
+
+ friend Ostream& operator<<
+ (
+ Ostream&,
+ const wallBoundedStreamLineParticle&
+ );
+
+ friend Ostream& operator<<
+ (
+ Ostream&,
+ const InfoProxy<wallBoundedStreamLineParticle>&
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
new file mode 100644
index 0000000000000000000000000000000000000000..dbe7ce827665f75efc529d49e258af19570ea7e6
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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 "wallBoundedStreamLineParticleCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTemplateTypeNameAndDebug(Cloud<wallBoundedStreamLineParticle>, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud
+(
+ const polyMesh& mesh,
+ const word& cloudName,
+ bool readFields
+)
+:
+ Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, false)
+{
+ if (readFields)
+ {
+ wallBoundedStreamLineParticle::readFields(*this);
+ }
+}
+
+
+Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud
+(
+ const polyMesh& mesh,
+ const word& cloudName,
+ const IDLList<wallBoundedStreamLineParticle>& particles
+)
+:
+ Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, particles)
+{}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H
new file mode 100644
index 0000000000000000000000000000000000000000..b91e9b5d8fa32007dc01cfe6e0426a9199027449
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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::wallBoundedStreamLineParticleCloud
+
+Description
+ A Cloud of streamLine particles
+
+SourceFiles
+ streamLineCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef streamLineCloud_H
+#define streamLineCloud_H
+
+#include "Cloud.H"
+#include "wallBoundedStreamLineParticle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class streamLineCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLineParticleCloud
+:
+ public Cloud<wallBoundedStreamLineParticle>
+{
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ wallBoundedStreamLineParticleCloud
+ (
+ const wallBoundedStreamLineParticleCloud&
+ );
+
+ //- Disallow default bitwise assignment
+ void operator=(const wallBoundedStreamLineParticleCloud&);
+
+
+public:
+
+ //- Type of parcel the cloud was instantiated for
+ typedef wallBoundedStreamLineParticle parcelType;
+
+ // Constructors
+
+ //- Construct given mesh
+ wallBoundedStreamLineParticleCloud
+ (
+ const polyMesh&,
+ const word& cloudName = "defaultCloud",
+ bool readFields = true
+ );
+
+ //- Construct from mesh, cloud name, and a list of particles
+ wallBoundedStreamLineParticleCloud
+ (
+ const polyMesh& mesh,
+ const word& cloudName,
+ const IDLList<wallBoundedStreamLineParticle>& particles
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //