ENH: nearWallFields: use tracking

This commit is contained in:
mattijs
2013-11-04 16:33:24 +00:00
parent c9386b4679
commit 1864c92db5
6 changed files with 832 additions and 63 deletions

View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "findCellParticle.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::findCellParticle::findCellParticle
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const point& end,
const label data
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
end_(end),
data_(data)
{}
Foam::findCellParticle::findCellParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
)
:
particle(mesh, is, readFields)
{
if (readFields)
{
if (is.format() == IOstream::ASCII)
{
is >> end_;
data_ = readLabel(is);
}
else
{
is.read
(
reinterpret_cast<char*>(&end_),
sizeof(end_) + sizeof(data_)
);
}
}
// Check state of Istream
is.check
(
"findCellParticle::findCellParticle"
"(const Cloud<findCellParticle>&, Istream&, bool)"
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::findCellParticle::move
(
trackingData& td,
const scalar maxTrackLen
)
{
td.switchProcessor = false;
td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*maxTrackLen;
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
{
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
dt *= trackToFace(end_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/maxTrackLen;
}
if (tEnd < SMALL || !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.
td.cellToData()[cell()].append(data());
td.cellToEnd()[cell()].append(position());
}
return td.keepParticle;
}
bool Foam::findCellParticle::hitPatch
(
const polyPatch&,
trackingData& td,
const label patchI,
const scalar trackFraction,
const tetIndices& tetIs
)
{
return false;
}
void Foam::findCellParticle::hitWedgePatch
(
const wedgePolyPatch&,
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::findCellParticle::hitSymmetryPatch
(
const symmetryPolyPatch&,
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::findCellParticle::hitCyclicPatch
(
const cyclicPolyPatch&,
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::findCellParticle::hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
)
{
// Remove particle
td.switchProcessor = true;
}
void Foam::findCellParticle::hitWallPatch
(
const wallPolyPatch& wpp,
trackingData& td,
const tetIndices&
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::findCellParticle::hitPatch
(
const polyPatch& wpp,
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
{
if (os.format() == IOstream::ASCII)
{
os << static_cast<const particle&>(p)
<< token::SPACE << p.end_
<< token::SPACE << p.data_;
}
else
{
os << static_cast<const particle&>(p);
os.write
(
reinterpret_cast<const char*>(&p.end_),
sizeof(p.end_) + sizeof(p.data_)
);
}
// Check state of Ostream
os.check("Ostream& operator<<(Ostream&, const findCellParticle&)");
return os;
}
// ************************************************************************* //

View File

@ -0,0 +1,259 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::findCellParticle
Description
Particle class that finds cells by tracking
SourceFiles
findCellParticle.C
\*---------------------------------------------------------------------------*/
#ifndef findCellParticle_H
#define findCellParticle_H
#include "particle.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class findCellParticleCloud;
/*---------------------------------------------------------------------------*\
Class findCellParticle Declaration
\*---------------------------------------------------------------------------*/
class findCellParticle
:
public particle
{
// Private data
//- end point to track to
point end_;
//- passive data
label data_;
public:
friend class Cloud<findCellParticle>;
//- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<Cloud<findCellParticle> >
{
labelListList& cellToData_;
List<List<point> >& cellToEnd_;
public:
// Constructors
trackingData
(
Cloud<findCellParticle>& cloud,
labelListList& cellToData,
List<List<point> >& cellToEnd
)
:
particle::TrackingData<Cloud<findCellParticle> >(cloud),
cellToData_(cellToData),
cellToEnd_(cellToEnd)
{}
// Member functions
labelListList& cellToData()
{
return cellToData_;
}
List<List<point> >& cellToEnd()
{
return cellToEnd_;
}
};
// Constructors
//- Construct from components
findCellParticle
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const point& end,
const label data
);
//- Construct from Istream
findCellParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields = true
);
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new findCellParticle(*this));
}
//- Factory class to read-construct particles used for
// parallel transfer
class iNew
{
const polyMesh& mesh_;
public:
iNew(const polyMesh& mesh)
:
mesh_(mesh)
{}
autoPtr<findCellParticle> operator()(Istream& is) const
{
return autoPtr<findCellParticle>
(
new findCellParticle(mesh_, is, true)
);
}
};
// Member Functions
//- point to track to
const point& end() const
{
return end_;
}
//- transported label
label data() const
{
return data_;
}
// Tracking
//- Track all particles to their end point
bool move(trackingData&, const scalar);
//- 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
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const findCellParticle&);
};
template<>
inline bool contiguous<findCellParticle>()
{
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,42 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "findCellParticle.H"
#include "Cloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTemplateTypeNameAndDebug(Cloud<findCellParticle>, 0);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -25,6 +25,9 @@ License
#include "nearWallFields.H" #include "nearWallFields.H"
#include "wordReList.H" #include "wordReList.H"
#include "findCellParticle.H"
#include "mappedPatchBase.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -34,6 +37,189 @@ defineTypeNameAndDebug(nearWallFields, 0);
} }
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::nearWallFields::calcAddressing()
{
const fvMesh& mesh = refCast<const fvMesh>(obr_);
// Count number of faces
label nPatchFaces = 0;
forAllConstIter(labelHashSet, patchSet_, iter)
{
label patchI = iter.key();
nPatchFaces += mesh.boundary()[patchI].size();
}
// Global indexing
globalIndex globalCells(mesh.nCells());
globalIndex globalWalls(nPatchFaces);
if (debug)
{
Info<< "nearWallFields::calcAddressing() :"
<< " nPatchFaces:" << globalWalls.size() << endl;
}
// Construct cloud
Cloud<findCellParticle> cloud(mesh, IDLList<findCellParticle>());
// Add particles to track to sample locations
nPatchFaces = 0;
forAllConstIter(labelHashSet, patchSet_, iter)
{
label patchI = iter.key();
const fvPatch& patch = mesh.boundary()[patchI];
vectorField nf = patch.nf();
vectorField faceCellCentres = patch.patch().faceCellCentres();
forAll(patch, patchFaceI)
{
label meshFaceI = patch.start()+patchFaceI;
// Find starting point on face (since faceCentre might not
// be on face-diagonal decomposition)
pointIndexHit startInfo
(
mappedPatchBase::facePoint
(
mesh,
meshFaceI,
polyMesh::FACEDIAGTETS
)
);
point start;
if (startInfo.hit())
{
start = startInfo.hitPoint();
}
else
{
// Fallback: start tracking from neighbouring cell centre
start = faceCellCentres[patchFaceI];
}
const point end = start-distance_*nf[patchFaceI];
// Find tet for starting location
label cellI = -1;
label tetFaceI = -1;
label tetPtI = -1;
mesh.findCellFacePt(start, cellI, tetFaceI, tetPtI);
// Add to cloud. Add originating face as passive data
cloud.addParticle
(
new findCellParticle
(
mesh,
start,
cellI,
tetFaceI,
tetPtI,
end,
globalWalls.toGlobal(nPatchFaces) // passive data
)
);
nPatchFaces++;
}
}
if (debug)
{
// Dump particles
OBJstream str
(
mesh.time().path()
/"wantedTracks_" + mesh.time().timeName() + ".obj"
);
Info<< "nearWallFields::calcAddressing() :"
<< "Dumping tracks to " << str.name() << endl;
forAllConstIter(Cloud<findCellParticle>, cloud, iter)
{
const findCellParticle& tp = iter();
str.write(linePointRef(tp.position(), tp.end()));
}
}
// Per cell: empty or global wall index and end location
cellToWalls_.setSize(mesh.nCells());
cellToSamples_.setSize(mesh.nCells());
// Database to pass into findCellParticle::move
findCellParticle::trackingData td(cloud, cellToWalls_, cellToSamples_);
// Track all particles to their end position.
scalar maxTrackLen = 2.0*mesh.bounds().mag();
//Debug: collect start points
pointField start;
if (debug)
{
start.setSize(nPatchFaces);
nPatchFaces = 0;
forAllConstIter(Cloud<findCellParticle>, cloud, iter)
{
const findCellParticle& tp = iter();
start[nPatchFaces++] = tp.position();
}
}
cloud.move(td, maxTrackLen);
// Rework cell-to-globalpatchface into a map
List<Map<label> > compactMap;
getPatchDataMapPtr_.reset
(
new mapDistribute
(
globalWalls,
cellToWalls_,
compactMap
)
);
// Debug: dump resulting tracks
if (debug)
{
getPatchDataMapPtr_().distribute(start);
{
OBJstream str
(
mesh.time().path()
/"obtainedTracks_" + mesh.time().timeName() + ".obj"
);
Info<< "nearWallFields::calcAddressing() :"
<< "Dumping obtained to " << str.name() << endl;
forAll(cellToWalls_, cellI)
{
const List<point>& ends = cellToSamples_[cellI];
const labelList& cData = cellToWalls_[cellI];
forAll(cData, i)
{
str.write(linePointRef(ends[i], start[cData[i]]));
}
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::nearWallFields::nearWallFields Foam::nearWallFields::nearWallFields
@ -127,16 +313,11 @@ void Foam::nearWallFields::read(const dictionary& dict)
reverseFieldMap_.insert(sampleFldName, fldName); reverseFieldMap_.insert(sampleFldName, fldName);
} }
Info<< type() << " " << name_ << ": Creating " << fieldMap_.size() Info<< type() << " " << name_ << ": Sampling " << fieldMap_.size()
<< " fields" << endl; << " fields" << endl;
createFields(vsf_); // Do analysis
createFields(vvf_); calcAddressing();
createFields(vSpheretf_);
createFields(vSymmtf_);
createFields(vtf_);
Info<< endl;
} }
} }
@ -147,15 +328,6 @@ void Foam::nearWallFields::execute()
{ {
Info<< "nearWallFields:execute()" << endl; Info<< "nearWallFields:execute()" << endl;
} }
//if (active_)
//{
// sampleFields(vsf_);
// sampleFields(vvf_);
// sampleFields(vSpheretf_);
// sampleFields(vSymmtf_);
// sampleFields(vtf_);
//}
} }
@ -165,8 +337,6 @@ void Foam::nearWallFields::end()
{ {
Info<< "nearWallFields:end()" << endl; Info<< "nearWallFields:end()" << endl;
} }
// Update fields
execute();
} }
@ -186,6 +356,28 @@ void Foam::nearWallFields::write()
// Do nothing // Do nothing
if (active_) if (active_)
{ {
if
(
fieldMap_.size()
&& vsf_.empty()
&& vvf_.empty()
&& vSpheretf_.empty()
&& vSymmtf_.empty()
&& vtf_.empty()
)
{
Info<< type() << " " << name_ << ": Creating " << fieldMap_.size()
<< " fields" << endl;
createFields(vsf_);
createFields(vvf_);
createFields(vSpheretf_);
createFields(vSymmtf_);
createFields(vtf_);
Info<< endl;
}
Info<< type() << " " << name_ << " output:" << nl; Info<< type() << " " << name_ << " output:" << nl;
Info<< " Writing sampled fields to " << obr_.time().timeName() Info<< " Writing sampled fields to " << obr_.time().timeName()

View File

@ -76,6 +76,7 @@ SourceFiles
#include "OFstream.H" #include "OFstream.H"
#include "volFields.H" #include "volFields.H"
#include "Tuple2.H" #include "Tuple2.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -122,21 +123,32 @@ protected:
//- From resulting back to original field //- From resulting back to original field
HashTable<word> reverseFieldMap_; HashTable<word> reverseFieldMap_;
//- Locally constructed fields
PtrList<volScalarField> vsf_; // Calculated addressing
PtrList<volVectorField> vvf_;
PtrList<volSphericalTensorField> vSpheretf_; //- From cell to seed patch faces
PtrList<volSymmTensorField> vSymmtf_; labelListList cellToWalls_;
PtrList<volTensorField> vtf_;
//- From cell to tracked end point
List<List<point> > cellToSamples_;
//- Map from cell based data back to patch based data
autoPtr<mapDistribute> getPatchDataMapPtr_;
// Locally constructed fields
PtrList<volScalarField> vsf_;
PtrList<volVectorField> vvf_;
PtrList<volSphericalTensorField> vSpheretf_;
PtrList<volSymmTensorField> vSymmtf_;
PtrList<volTensorField> vtf_;
// Protected Member Functions // Protected Member Functions
//- Disallow default bitwise copy construct //- Calculate addressing from cells back to patch faces
nearWallFields(const nearWallFields&); void calcAddressing();
//- Disallow default bitwise assignment
void operator=(const nearWallFields&);
template<class Type> template<class Type>
void createFields void createFields
@ -144,6 +156,14 @@ protected:
PtrList<GeometricField<Type, fvPatchField, volMesh> >& PtrList<GeometricField<Type, fvPatchField, volMesh> >&
) const; ) const;
//- Override boundary fields with sampled values
template<class Type>
void sampleBoundaryField
(
const interpolationCellPoint<Type>& interpolator,
GeometricField<Type, fvPatchField, volMesh>& fld
) const;
template<class Type> template<class Type>
void sampleFields void sampleFields
( (
@ -151,6 +171,14 @@ protected:
) const; ) const;
private:
//- Disallow default bitwise copy construct
nearWallFields(const nearWallFields&);
//- Disallow default bitwise assignment
void operator=(const nearWallFields&);
public: public:
//- Runtime type information //- Runtime type information

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -24,8 +24,6 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "nearWallFields.H" #include "nearWallFields.H"
#include "mappedFieldFvPatchFields.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,35 +62,8 @@ void Foam::nearWallFields::createFields
io.rename(sampleFldName); io.rename(sampleFldName);
sflds.set(sz, new vfType(io, fld)); sflds.set(sz, new vfType(io, fld));
vfType& sampleFld = sflds[sz];
// Reset the bcs to be mapped Info<< " created " << sflds[sz].name() << " to sample "
forAllConstIter(labelHashSet, patchSet_, iter)
{
label patchI = iter.key();
sampleFld.boundaryField().set
(
patchI,
new mappedFieldFvPatchField<Type>
(
sampleFld.mesh().boundary()[patchI],
sampleFld.dimensionedInternalField(),
sampleFld.mesh().name(),
mappedPatchBase::NEARESTCELL,
word::null, // samplePatch
-distance_,
sampleFld.name(), // fieldName
false, // setAverage
pTraits<Type>::zero, // average
interpolationCellPoint<Type>::typeName
)
);
}
Info<< " created " << sampleFld.name() << " to sample "
<< fld.name() << endl; << fld.name() << endl;
} }
} }
@ -100,6 +71,53 @@ void Foam::nearWallFields::createFields
} }
template<class Type>
void Foam::nearWallFields::sampleBoundaryField
(
const interpolationCellPoint<Type>& interpolator,
GeometricField<Type, fvPatchField, volMesh>& fld
) const
{
// Construct flat fields for all patch faces to be sampled
Field<Type> sampledValues(getPatchDataMapPtr_().constructSize());
forAll(cellToWalls_, cellI)
{
const labelList& cData = cellToWalls_[cellI];
forAll(cData, i)
{
const point& samplePt = cellToSamples_[cellI][i];
sampledValues[cData[i]] = interpolator.interpolate(samplePt, cellI);
}
}
// Send back sampled values to patch faces
getPatchDataMapPtr_().reverseDistribute
(
getPatchDataMapPtr_().constructSize(),
sampledValues
);
// Pick up data
label nPatchFaces = 0;
forAllConstIter(labelHashSet, patchSet_, iter)
{
label patchI = iter.key();
fvPatchField<Type>& pfld = fld.boundaryField()[patchI];
Field<Type> newFld(pfld.size());
forAll(pfld, i)
{
newFld[i] = sampledValues[nPatchFaces++];
}
pfld == newFld;
}
}
template<class Type> template<class Type>
void Foam::nearWallFields::sampleFields void Foam::nearWallFields::sampleFields
( (
@ -115,8 +133,12 @@ void Foam::nearWallFields::sampleFields
// Take over internal and boundary values // Take over internal and boundary values
sflds[i] == fld; sflds[i] == fld;
// Evaluate to update the mapped
sflds[i].correctBoundaryConditions(); // Construct interpolation method
interpolationCellPoint<Type> interpolator(fld);
// Override sampled values
sampleBoundaryField(interpolator, sflds[i]);
} }
} }