zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'particleInteractions' of ssh://hunt//home/noisy3/OpenFOAM/OpenFOAM-dev into particleInteractions

Browse Source
This commit is contained in:
graham
2010-03-16 22:16:54 +00:00
parent 33bc146b88 5b7ef77afa
commit 8d5ecca7ae
82 changed files with 4889 additions and 1109 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C
View File

@ -40,7 +40,7 @@ template <class Type>
const Foam::scalar Foam::FaceCellWave<Type>::geomTol_ = 1e-6;
template <class Type>
const Foam::scalar Foam::FaceCellWave<Type>::propagationTol_ = 0.01;
Foam::scalar Foam::FaceCellWave<Type>::propagationTol_ = 0.01;
// Write to ostream
template <class Type>

2
src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H
View File

@ -252,7 +252,7 @@ class FaceCellWave
// Private static data
static const scalar geomTol_;
static const scalar propagationTol_;
static scalar propagationTol_;
public:

12
src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C
View File

@ -38,7 +38,17 @@ Foam::IPstream::IPstream
)
:
Pstream(commsType, bufSize),
UIPstream(commsType, fromProcNo, buf_, externalBufPosition_),
UIPstream
(
commsType,
fromProcNo,
buf_,
externalBufPosition_,
UPstream::msgType(), // tag
false, // do not clear buf_ if at end
format,
version
),
externalBufPosition_(0)
{}

15
src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C
View File

@ -77,6 +77,21 @@ inline void Foam::UIPstream::readFromBuffer
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::UIPstream::~UIPstream()
{
if (clearAtEnd_ && eof())
{
if (debug)
{
Pout<< "UIPstream::~UIPstream() : clearing externalBuf_ of size "
<< externalBuf_.size() << endl;
}
externalBuf_.clearStorage();
}
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::Istream& Foam::UIPstream::read(token& t)

8
src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H
View File

@ -67,6 +67,8 @@ class UIPstream
const int tag_;
const bool clearAtEnd_;
label messageSize_;
@ -96,6 +98,7 @@ public:
DynamicList<char>& externalBuf,
label& externalBufPosition,
const int tag = UPstream::msgType(),
const bool clearAtEnd = false, // destroy externalBuf if at end
streamFormat format=BINARY,
versionNumber version=currentVersion
);
@ -104,6 +107,11 @@ public:
UIPstream(const int fromProcNo, PstreamBuffers&);
// Destructor
~UIPstream();
// Member functions
// Inquiry

16
src/OpenFOAM/db/Time/Time.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -625,21 +625,25 @@ void Foam::Time::setDeltaT(const scalar deltaT)
Foam::TimeState Foam::Time::subCycle(const label nSubCycles)
{
subCycling_ = true;
prevTimeState_.set(new TimeState(*this));
TimeState ts = *this;
setTime(*this - deltaT(), (timeIndex() - 1)*nSubCycles);
deltaT_ /= nSubCycles;
deltaT0_ /= nSubCycles;
deltaTSave_ = deltaT0_;
return ts;
return prevTimeState();
}
void Foam::Time::endSubCycle(const TimeState& ts)
void Foam::Time::endSubCycle()
{
subCycling_ = false;
TimeState::operator=(ts);
if (subCycling_)
{
subCycling_ = false;
TimeState::operator=(prevTimeState());
prevTimeState_.clear();
}
}

21
src/OpenFOAM/db/Time/Time.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,6 +127,9 @@ protected:
//- Is the time currently being sub-cycled?
bool subCycling_;
//- If time is being sub-cycled this is the previous TimeState
autoPtr<TimeState> prevTimeState_;
//- Time directory name format
static fmtflags format_;
@ -347,6 +350,18 @@ public:
return functionObjects_;
}
//- Return true if time currently being sub-cycled, otherwise false
bool subCycling() const
{
return subCycling_;
}
//- Return previous TimeState if time is being sub-cycled
const TimeState& prevTimeState() const
{
return prevTimeState_();
}
// Check
@ -427,8 +442,8 @@ public:
//- Set time to sub-cycle for the given number of steps
virtual TimeState subCycle(const label nSubCycles);
//- Reset time after sub-cycling back to given TimeState
virtual void endSubCycle(const TimeState&);
//- Reset time after sub-cycling back to previous TimeState
virtual void endSubCycle();
//- Return non-const access to the list of function objects
functionObjectList& functionObjects()

9
src/OpenFOAM/db/Time/subCycleTime.C
View File

@ -32,9 +32,10 @@ Foam::subCycleTime::subCycleTime(Time& t, const label nSubCycles)
:
time_(t),
nSubCycles_(nSubCycles),
subCycleIndex_(0),
initialTimeState_(time_.subCycle(nSubCycles_))
{}
subCycleIndex_(0)
{
time_.subCycle(nSubCycles_);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
@ -55,7 +56,7 @@ bool Foam::subCycleTime::end() const
void Foam::subCycleTime::endSubCycle()
{
time_.endSubCycle(initialTimeState_);
time_.endSubCycle();
}

1
src/OpenFOAM/db/Time/subCycleTime.H
View File

@ -55,7 +55,6 @@ class subCycleTime
label nSubCycles_;
label subCycleIndex_;
TimeState initialTimeState_;
public:

3
src/Pstream/dummy/UIPread.C
View File

@ -38,6 +38,7 @@ Foam::UIPstream::UIPstream
DynamicList<char>& externalBuf,
label& externalBufPosition,
const int tag,
const bool clearAtEnd,
streamFormat format,
versionNumber version
)
@ -48,6 +49,7 @@ Foam::UIPstream::UIPstream
externalBuf_(externalBuf),
externalBufPosition_(externalBufPosition),
tag_(tag),
clearAtEnd_(clearAtEnd),
messageSize_(0)
{
notImplemented
@ -59,6 +61,7 @@ Foam::UIPstream::UIPstream
"DynamicList<char>&,"
"label&,"
"const int tag,"
"const bool,"
"streamFormat, versionNumber"
")"
);

3
src/Pstream/mpi/UIPread.C
View File

@ -42,6 +42,7 @@ Foam::UIPstream::UIPstream
DynamicList<char>& externalBuf,
label& externalBufPosition,
const int tag,
const bool clearAtEnd,
streamFormat format,
versionNumber version
)
@ -52,6 +53,7 @@ Foam::UIPstream::UIPstream
externalBuf_(externalBuf),
externalBufPosition_(externalBufPosition),
tag_(tag),
clearAtEnd_(clearAtEnd),
messageSize_(0)
{
setOpened();
@ -125,6 +127,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
externalBuf_(buffers.recvBuf_[fromProcNo]),
externalBufPosition_(buffers.recvBufPos_[fromProcNo]),
tag_(buffers.tag_),
clearAtEnd_(true),
messageSize_(0)
{
if (commsType() != UPstream::scheduled && !buffers.finishedSendsCalled_)

1
src/finiteVolume/Make/files
View File

@ -242,7 +242,6 @@ $(limitedSchemes)/SuperBee/SuperBee.C
$(limitedSchemes)/QUICK/QUICK.C
$(limitedSchemes)/MUSCL/MUSCL.C
$(limitedSchemes)/UMIST/UMIST.C
$(limitedSchemes)/MC/MC.C
$(limitedSchemes)/Phi/Phi.C
$(limitedSchemes)/filteredLinear/filteredLinear.C
$(limitedSchemes)/filteredLinear2/filteredLinear2.C

14
src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,6 +34,7 @@ License
#include "fvcSurfaceIntegrate.H"
#include "slicedSurfaceFields.H"
#include "syncTools.H"
#include "fvm.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -107,7 +108,7 @@ void Foam::MULES::explicitSolve
{
psiIf =
(
mesh.V0()*rho.oldTime()*psi0/(deltaT*mesh.V())
mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
+ Su.field()
- psiIf
)/(rho/deltaT - Sp.field());
@ -325,7 +326,8 @@ void Foam::MULES::limiter
const unallocLabelList& owner = mesh.owner();
const unallocLabelList& neighb = mesh.neighbour();
const scalarField& V = mesh.V();
tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
const scalar deltaT = mesh.time().deltaTValue();
const scalarField& phiBDIf = phiBD;
@ -452,14 +454,16 @@ void Foam::MULES::limiter
if (mesh.moving())
{
tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
psiMaxn =
V*((rho/deltaT - Sp)*psiMaxn - Su)
- (mesh.V0()/deltaT)*rho.oldTime()*psi0
- (V0()/deltaT)*rho.oldTime()*psi0
+ sumPhiBD;
psiMinn =
V*(Su - (rho/deltaT - Sp)*psiMinn)
+ (mesh.V0()/deltaT)*rho.oldTime()*psi0
+ (V0/deltaT)*rho.oldTime()*psi0
- sumPhiBD;
}
else

6
src/finiteVolume/fvMesh/fvMesh.H
View File

@ -269,6 +269,12 @@ public:
//- Return old-old-time cell volumes
const DimensionedField<scalar, volMesh>& V00() const;
//- Return sub-cycle cell volumes
tmp<DimensionedField<scalar, volMesh> > Vsc() const;
//- Return sub-cycl old-time cell volumes
tmp<DimensionedField<scalar, volMesh> > Vsc0() const;
//- Return cell face area vectors
const surfaceVectorField& Sf() const;

57
src/finiteVolume/fvMesh/fvMeshGeometry.C
View File

@ -285,6 +285,63 @@ const volScalarField::DimensionedInternalField& fvMesh::V00() const
}
tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
{
if (moving() && time().subCycling())
{
const TimeState& ts = time();
const TimeState& ts0 = time().prevTimeState();
scalar tFrac =
(
ts.value() - (ts0.value() - ts0.deltaTValue())
)/ts0.deltaTValue();
if (tFrac < (1 - SMALL))
{
return V0() + tFrac*(V() - V0());
}
else
{
return V();
}
}
else
{
return V();
}
}
tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
{
if (moving() && time().subCycling())
{
const TimeState& ts = time();
const TimeState& ts0 = time().prevTimeState();
scalar t0Frac =
(
(ts.value() - ts.deltaTValue())
- (ts0.value() - ts0.deltaTValue())
)/ts0.deltaTValue();
if (t0Frac > SMALL)
{
return V0() + t0Frac*(V() - V0());
}
else
{
return V0();
}
}
else
{
return V0();
}
}
const surfaceVectorField& fvMesh::Sf() const
{
if (!SfPtr_)

94
src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/MC/MC.H
View File

@ -1,94 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::MCLimiter
Description
Class with limiter function which returns the limiter for the
monotonised centred differencing scheme based on r obtained from
the LimiterFunc class.
Used in conjunction with the template class LimitedScheme.
SourceFiles
MC.C
\*---------------------------------------------------------------------------*/
#ifndef MC_H
#define MC_H
#include "vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class MCLimiter Declaration
\*---------------------------------------------------------------------------*/
template<class LimiterFunc>
class MCLimiter
:
public LimiterFunc
{
public:
MCLimiter(Istream&)
{}
scalar limiter
(
const scalar cdWeight,
const scalar faceFlux,
const typename LimiterFunc::phiType& phiP,
const typename LimiterFunc::phiType& phiN,
const typename LimiterFunc::gradPhiType& gradcP,
const typename LimiterFunc::gradPhiType& gradcN,
const vector& d
) const
{
scalar r = LimiterFunc::r
(
faceFlux, phiP, phiN, gradcP, gradcN, d
);
return max(min(min(2*r, 0.5*(1 + r)), 2), 0);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/fvMotionSolver/Make/files
View File

@ -1,6 +1,8 @@
fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C
fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

567
src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C Normal file
View File

@ -0,0 +1,567 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "displacementLayeredMotionFvMotionSolver.H"
#include "addToRunTimeSelectionTable.H"
#include "pointEdgeStructuredWalk.H"
#include "pointFields.H"
#include "PointEdgeWave.H"
#include "syncTools.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(displacementLayeredMotionFvMotionSolver, 0);
addToRunTimeSelectionTable
(
fvMotionSolver,
displacementLayeredMotionFvMotionSolver,
dictionary
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::displacementLayeredMotionFvMotionSolver::calcZoneMask
(
const label cellZoneI,
PackedBoolList& isZonePoint,
PackedBoolList& isZoneEdge
) const
{
if (cellZoneI == -1)
{
isZonePoint.setSize(mesh().nPoints());
isZonePoint = 1;
isZoneEdge.setSize(mesh().nEdges());
isZoneEdge = 1;
}
else
{
const cellZone& cz = mesh().cellZones()[cellZoneI];
label nPoints = 0;
forAll(cz, i)
{
const labelList& cPoints = mesh().cellPoints(cz[i]);
forAll(cPoints, cPointI)
{
if (!isZonePoint[cPoints[cPointI]])
{
isZonePoint[cPoints[cPointI]] = 1;
nPoints++;
}
}
}
syncTools::syncPointList
(
mesh(),
isZonePoint,
orEqOp<unsigned int>(),
0
);
// Mark edge inside cellZone
label nEdges = 0;
forAll(cz, i)
{
const labelList& cEdges = mesh().cellEdges(cz[i]);
forAll(cEdges, cEdgeI)
{
if (!isZoneEdge[cEdges[cEdgeI]])
{
isZoneEdge[cEdges[cEdgeI]] = 1;
nEdges++;
}
}
}
syncTools::syncEdgeList
(
mesh(),
isZoneEdge,
orEqOp<unsigned int>(),
0
);
Info<< "On cellZone " << cz.name()
<< " marked " << returnReduce(nPoints, sumOp<label>())
<< " points and " << returnReduce(nEdges, sumOp<label>())
<< " edges." << endl;
}
}
// Find distance to starting point
void Foam::displacementLayeredMotionFvMotionSolver::walkStructured
(
const label cellZoneI,
const PackedBoolList& isZonePoint,
const PackedBoolList& isZoneEdge,
const labelList& seedPoints,
const vectorField& seedData,
scalarField& distance,
vectorField& data
) const
{
const pointField& points = mesh().points();
List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
forAll(seedPoints, i)
{
seedInfo[i] = pointEdgeStructuredWalk
(
true,
points[seedPoints[i]],
0.0,
seedData[i]
);
}
// Current info on points
List<pointEdgeStructuredWalk> allPointInfo(mesh().nPoints());
// Mark points inside cellZone
forAll(isZonePoint, pointI)
{
if (isZonePoint[pointI])
{
allPointInfo[pointI].inZone() = true;
}
}
// Current info on edges
List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
// Mark edges inside cellZone
forAll(isZoneEdge, edgeI)
{
if (isZoneEdge[edgeI])
{
allEdgeInfo[edgeI].inZone() = true;
}
}
// Walk
PointEdgeWave<pointEdgeStructuredWalk> wallCalc
(
mesh(),
seedPoints,
seedInfo,
allPointInfo,
allEdgeInfo,
mesh().globalData().nTotalPoints() // max iterations
);
// Extract distance and passive data
forAll(allPointInfo, pointI)
{
if (isZonePoint[pointI])
{
distance[pointI] = allPointInfo[pointI].dist();
data[pointI] = allPointInfo[pointI].data();
}
}
}
// Evaluate faceZone patch
Foam::tmp<Foam::vectorField>
Foam::displacementLayeredMotionFvMotionSolver::faceZoneEvaluate
(
const faceZone& fz,
const labelList& meshPoints,
const dictionary& dict,
const PtrList<pointVectorField>& patchDisp,
const label patchI
) const
{
tmp<vectorField> tfld(new vectorField(meshPoints.size()));
vectorField& fld = tfld();
const word& type = dict.lookup("type");
if (type == "fixedValue")
{
fld = vectorField("value", dict, meshPoints.size());
}
else if (type == "timeVaryingUniformFixedValue")
{
interpolationTable<vector> timeSeries(dict);
fld = timeSeries(mesh().time().timeOutputValue());
}
else if (type == "slip")
{
if ((patchI % 2) != 1)
{
FatalIOErrorIn
(
"displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
*this
) << "slip can only be used on second faceZonePatch of pair."
<< "FaceZone:" << fz.name()
<< exit(FatalIOError);
}
// Use field set by previous bc
fld = vectorField(patchDisp[patchI-1], meshPoints);
}
else if (type == "follow")
{
// Only on boundary faces - follow boundary conditions
fld = vectorField(pointDisplacement_, meshPoints);
}
else
{
FatalIOErrorIn
(
"displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
*this
) << "Unknown faceZonePatch type " << type << " for faceZone "
<< fz.name() << exit(FatalIOError);
}
return tfld;
}
void Foam::displacementLayeredMotionFvMotionSolver::cellZoneSolve
(
const label cellZoneI,
const dictionary& zoneDict
)
{
PackedBoolList isZonePoint(mesh().nPoints());
PackedBoolList isZoneEdge(mesh().nEdges());
calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
const dictionary& patchesDict = zoneDict.subDict("boundaryField");
if (patchesDict.size() != 2)
{
FatalIOErrorIn
(
"displacementLayeredMotionFvMotionSolver::"
"correctBoundaryConditions(..)",
*this
) << "Can only handle 2 faceZones (= patches) per cellZone. "
<< " cellZone:" << cellZoneI
<< " patches:" << patchesDict.toc()
<< exit(FatalIOError);
}
PtrList<scalarField> patchDist(patchesDict.size());
PtrList<pointVectorField> patchDisp(patchesDict.size());
// Allocate the fields
label patchI = 0;
forAllConstIter(dictionary, patchesDict, patchIter)
{
const word& faceZoneName = patchIter().keyword();
label zoneI = mesh().faceZones().findZoneID(faceZoneName);
if (zoneI == -1)
{
FatalIOErrorIn
(
"displacementLayeredMotionFvMotionSolver::"
"correctBoundaryConditions(..)",
*this
) << "Cannot find faceZone " << faceZoneName
<< endl << "Valid zones are " << mesh().faceZones().names()
<< exit(FatalIOError);
}
// Determine the points of the faceZone within the cellZone
const faceZone& fz = mesh().faceZones()[zoneI];
patchDist.set(patchI, new scalarField(mesh().nPoints()));
patchDisp.set
(
patchI,
new pointVectorField
(
IOobject
(
mesh().cellZones()[cellZoneI].name() + "_" + fz.name(),
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
pointDisplacement_ // to inherit the boundary conditions
)
);
patchI++;
}
// 'correctBoundaryConditions'
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Loops over all the faceZones and walks their boundary values
// Make sure we can pick up bc values from field
pointDisplacement_.correctBoundaryConditions();
patchI = 0;
forAllConstIter(dictionary, patchesDict, patchIter)
{
const word& faceZoneName = patchIter().keyword();
const dictionary& faceZoneDict = patchIter().dict();
// Determine the points of the faceZone within the cellZone
label zoneI = mesh().faceZones().findZoneID(faceZoneName);
const faceZone& fz = mesh().faceZones()[zoneI];
const labelList& fzMeshPoints = fz().meshPoints();
DynamicList<label> meshPoints(fzMeshPoints.size());
forAll(fzMeshPoints, i)
{
if (isZonePoint[fzMeshPoints[i]])
{
meshPoints.append(fzMeshPoints[i]);
}
}
// Get initial value for all the faceZone points
tmp<vectorField> tseed = faceZoneEvaluate
(
fz,
meshPoints,
faceZoneDict,
patchDisp,
patchI
);
Info<< "For cellZone:" << cellZoneI
<< " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
<< " have patchField:"
<< " max:" << gMax(tseed())
<< " min:" << gMin(tseed())
<< " avg:" << gAverage(tseed())
<< endl;
// Set distance and transported value
walkStructured
(
cellZoneI,
isZonePoint,
isZoneEdge,
meshPoints,
tseed,
patchDist[patchI],
patchDisp[patchI]
);
// Implement real bc.
patchDisp[patchI].correctBoundaryConditions();
//Info<< "Writing displacement for faceZone " << fz.name()
// << " to " << patchDisp[patchI].name() << endl;
//patchDisp[patchI].write();
// // Copy into pointDisplacement for other fields to use
// forAll(isZonePoint, pointI)
// {
// if (isZonePoint[pointI])
// {
// pointDisplacement_[pointI] = patchDisp[patchI][pointI];
// }
// }
// pointDisplacement_.correctBoundaryConditions();
patchI++;
}
// Solve
// ~~~~~
// solving the interior is just interpolating
// // Get normalised distance
// pointScalarField distance
// (
// IOobject
// (
// "distance",
// mesh().time().timeName(),
// mesh(),
// IOobject::NO_READ,
// IOobject::NO_WRITE,
// false
// ),
// pointMesh::New(mesh()),
// dimensionedScalar("distance", dimLength, 0.0)
// );
// forAll(distance, pointI)
// {
// if (isZonePoint[pointI])
// {
// scalar d1 = patchDist[0][pointI];
// scalar d2 = patchDist[1][pointI];
// if (d1+d2 > SMALL)
// {
// scalar s = d1/(d1+d2);
// distance[pointI] = s;
// }
// }
// }
// Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
// << endl;
// distance.write();
// Average
forAll(pointDisplacement_, pointI)
{
if (isZonePoint[pointI])
{
scalar d1 = patchDist[0][pointI];
scalar d2 = patchDist[1][pointI];
scalar s = d1/(d1+d2+VSMALL);
pointDisplacement_[pointI] =
(1-s)*patchDisp[0][pointI]
+ s*patchDisp[1][pointI];
}
}
pointDisplacement_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::displacementLayeredMotionFvMotionSolver::
displacementLayeredMotionFvMotionSolver
(
const polyMesh& mesh,
Istream& is
)
:
displacementFvMotionSolver(mesh, is),
pointDisplacement_
(
IOobject
(
"pointDisplacement",
fvMesh_.time().timeName(),
fvMesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
pointMesh::New(fvMesh_)
)
{
pointDisplacement_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::displacementLayeredMotionFvMotionSolver::
~displacementLayeredMotionFvMotionSolver()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::pointField>
Foam::displacementLayeredMotionFvMotionSolver::curPoints() const
{
tmp<pointField> tcurPoints
(
points0() + pointDisplacement_.internalField()
);
twoDCorrectPoints(tcurPoints());
// const pointField& pts = tcurPoints();
// forAll(pts, pointI)
// {
// Info<< " from:" << mesh().points()[pointI]
// << " to:" << pts[pointI]
// << endl;
// }
return tcurPoints;
}
void Foam::displacementLayeredMotionFvMotionSolver::solve()
{
const dictionary& ms = mesh().lookupObject<motionSolver>("dynamicMeshDict");
const dictionary& solverDict = ms.subDict(typeName + "Coeffs");
// Apply all regions (=cellZones)
const dictionary& regionDicts = solverDict.subDict("regions");
forAllConstIter(dictionary, regionDicts, regionIter)
{
const word& cellZoneName = regionIter().keyword();
const dictionary& regionDict = regionIter().dict();
label zoneI = mesh().cellZones().findZoneID(cellZoneName);
Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
<< endl;
if (zoneI == -1)
{
FatalIOErrorIn
(
"displacementLayeredMotionFvMotionSolver::solve(..)",
*this
) << "Cannot find cellZone " << cellZoneName
<< endl << "Valid zones are " << mesh().cellZones().names()
<< exit(FatalIOError);
}
cellZoneSolve(zoneI, regionDict);
}
}
void Foam::displacementLayeredMotionFvMotionSolver::updateMesh
(
const mapPolyMesh& mpm
)
{
displacementFvMotionSolver::updateMesh(mpm);
}
// ************************************************************************* //

187
src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H Normal file
View File

@ -0,0 +1,187 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::displacementLayeredMotionFvMotionSolver
Description
Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
structure of the mesh blocks and boundary conditions on these blocks.
Note: should not be an fvMotionSolver but just a motionSolver. Only here
so we can reuse displacementFvMotionSolver functionality (e.g. surface
following boundary conditions)
The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
per region (=cellZone) the boundary conditions on two opposing patches
(=faceZones). It then interpolates the boundary values using topological
walking so requires the cellZone to be a layered mesh.
The boundary conditions on faceZones are currently:
follow: the faceZone follows the overall mesh displacement.
Use this for faceZones on boundary faces (so it uses the
proper boundary conditions on the pointDisplacement).
fixedValue: fixed value.
timeVaryingUniformFixedValue: table-driven fixed value.
slip: the second of a pair of faceZones follows the tangential movement
specified by the first faceZone. (= removes the normal component).
SourceFiles
displacementLayeredMotionFvMotionSolver.C
\*---------------------------------------------------------------------------*/
#ifndef displacementLayeredMotionFvMotionSolver_H
#define displacementLayeredMotionFvMotionSolver_H
#include "displacementFvMotionSolver.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward class declarations
/*---------------------------------------------------------------------------*\
Class displacementLayeredMotionFvMotionSolver Declaration
\*---------------------------------------------------------------------------*/
class displacementLayeredMotionFvMotionSolver
:
public displacementFvMotionSolver
{
// Private data
//- Point motion field
mutable pointVectorField pointDisplacement_;
// Private Member Functions
void calcZoneMask
(
const label cellZoneI,
PackedBoolList& isZonePoint,
PackedBoolList& isZoneEdge
) const;
void walkStructured
(
const label cellZoneI,
const PackedBoolList& isZonePoint,
const PackedBoolList& isZoneEdge,
const labelList& seedPoints,
const vectorField& seedData,
scalarField& distance,
vectorField& data
) const;
tmp<vectorField> faceZoneEvaluate
(
const faceZone& fz,
const labelList& meshPoints,
const dictionary& dict,
const PtrList<pointVectorField>& patchDisp,
const label patchI
) const;
void cellZoneSolve
(
const label cellZoneI,
const dictionary& zoneDict
);
//- Disallow default bitwise copy construct
displacementLayeredMotionFvMotionSolver
(
const displacementLayeredMotionFvMotionSolver&
);
//- Disallow default bitwise assignment
void operator=(const displacementLayeredMotionFvMotionSolver&);
public:
//- Runtime type information
TypeName("displacementLayeredMotion");
// Constructors
//- Construct from polyMesh and data stream
displacementLayeredMotionFvMotionSolver
(
const polyMesh&,
Istream& msDataUnused
);
// Destructor
~displacementLayeredMotionFvMotionSolver();
// Member Functions
//- Return reference to the point motion displacement field
pointVectorField& pointDisplacement()
{
return pointDisplacement_;
}
//- Return const reference to the point motion displacement field
const pointVectorField& pointDisplacement() const
{
return pointDisplacement_;
}
//- Return point location obtained from the current motion field
virtual tmp<pointField> curPoints() const;
//- Solve for motion
virtual void solve();
//- Update topology
virtual void updateMesh(const mapPolyMesh&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

28
src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/MC/MC.C → src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
View File

@ -24,15 +24,31 @@ License
\*---------------------------------------------------------------------------*/
#include "LimitedScheme.H"
#include "MC.H"
#include "pointEdgeStructuredWalk.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
namespace Foam
Foam::Ostream& Foam::operator<<
(
Foam::Ostream& os,
const Foam::pointEdgeStructuredWalk& wDist
)
{
makeLimitedSurfaceInterpolationScheme(MC, MCLimiter)
makeLimitedVSurfaceInterpolationScheme(MCV, MCLimiter)
return os
<< wDist.inZone_ << wDist.previousPoint_
<< wDist.dist_ << wDist.data_;
}
Foam::Istream& Foam::operator>>
(
Foam::Istream& is,
Foam::pointEdgeStructuredWalk& wDist
)
{
return is
>> wDist.inZone_ >> wDist.previousPoint_
>> wDist.dist_ >> wDist.data_;
}
// ************************************************************************* //

224
src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H Normal file
View File

@ -0,0 +1,224 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::pointEdgeStructuredWalk
Description
Determines length of string of edges walked to point.
SourceFiles
pointEdgeStructuredWalkI.H
pointEdgeStructuredWalk.C
\*---------------------------------------------------------------------------*/
#ifndef pointEdgeStructuredWalk_H
#define pointEdgeStructuredWalk_H
#include "point.H"
#include "label.H"
#include "scalar.H"
#include "tensor.H"
#include "pTraits.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class polyPatch;
class polyMesh;
/*---------------------------------------------------------------------------*\
Class pointEdgeStructuredWalk Declaration
\*---------------------------------------------------------------------------*/
class pointEdgeStructuredWalk
{
// Private data
//- Is this point/edge inside the walk zone
bool inZone_;
//- Previuos point
point previousPoint_;
//- Sum of distance
scalar dist_;
//- Passive data
vector data_;
// Private Member Functions
//- Evaluate distance to point.
inline bool update
(
const point&,
const pointEdgeStructuredWalk& w2,
const scalar tol
);
//- Combine current with w2. Update distSqr, origin if w2 has smaller
// quantities and returns true.
inline bool update
(
const pointEdgeStructuredWalk& w2,
const scalar tol
);
public:
// Constructors
//- Construct null
inline pointEdgeStructuredWalk();
//- Construct from components
inline pointEdgeStructuredWalk
(
const bool,
const point&,
const scalar,
const vector&
);
// Member Functions
// Access
inline bool inZone() const;
inline bool& inZone();
inline const point& previousPoint() const;
inline scalar dist() const;
inline const vector& data() const;
// Needed by meshWave
//- Check whether still contains original (invalid) value.
inline bool valid() const;
//- Check for identical geometrical data. Used for cyclics checking.
inline bool sameGeometry
(
const pointEdgeStructuredWalk&,
const scalar tol
) const;
//- Convert origin to relative vector to leaving point
// (= point coordinate)
inline void leaveDomain
(
const polyPatch& patch,
const label patchPointI,
const point& pos
);
//- Convert relative origin to absolute by adding entering point
inline void enterDomain
(
const polyPatch& patch,
const label patchPointI,
const point& pos
);
//- Apply rotation matrix to origin
inline void transform(const tensor& rotTensor);
//- Influence of edge on point
inline bool updatePoint
(
const polyMesh& mesh,
const label pointI,
const label edgeI,
const pointEdgeStructuredWalk& edgeInfo,
const scalar tol
);
//- Influence of different value on same point.
// Merge new and old info.
inline bool updatePoint
(
const polyMesh& mesh,
const label pointI,
const pointEdgeStructuredWalk& newPointInfo,
const scalar tol
);
//- Influence of different value on same point.
// No information about current position whatsoever.
inline bool updatePoint
(
const pointEdgeStructuredWalk& newPointInfo,
const scalar tol
);
//- Influence of point on edge.
inline bool updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointI,
const pointEdgeStructuredWalk& pointInfo,
const scalar tol
);
// Member Operators
//Note: Used to determine whether to call update.
inline bool operator==(const pointEdgeStructuredWalk&) const;
inline bool operator!=(const pointEdgeStructuredWalk&) const;
// IOstream Operators
friend Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
friend Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "pointEdgeStructuredWalkI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

302
src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H Normal file
View File

@ -0,0 +1,302 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "polyMesh.H"
#include "transform.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Update this with w2.
inline bool Foam::pointEdgeStructuredWalk::update
(
const point& pt,
const pointEdgeStructuredWalk& w2,
const scalar tol
)
{
if (!valid())
{
// current not yet set. Walked from w2 to here (=pt)
dist_ = w2.dist_ + mag(pt-w2.previousPoint_);
previousPoint_ = pt;
data_ = w2.data_;
return true;
}
else
{
return false;
}
}
// Update this with w2 (information on same point)
inline bool Foam::pointEdgeStructuredWalk::update
(
const pointEdgeStructuredWalk& w2,
const scalar tol
)
{
if (!valid())
{
// current not yet set so use neighbour
dist_ = w2.dist_;
previousPoint_ = w2.previousPoint_;
data_ = w2.data_;
return true;
}
else
{
// already nearer to pt
return false;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Null constructor
inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
:
inZone_(false),
previousPoint_(vector::max),
dist_(0),
data_(vector::zero)
{}
// Construct from origin, distance
inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
(
const bool inZone,
const point& previousPoint,
const scalar dist,
const vector& data
)
:
inZone_(inZone),
previousPoint_(previousPoint),
dist_(dist),
data_(data)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline bool Foam::pointEdgeStructuredWalk::inZone() const
{
return inZone_;
}
inline bool& Foam::pointEdgeStructuredWalk::inZone()
{
return inZone_;
}
inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
{
return previousPoint_;
}
inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
{
return dist_;
}
inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
{
return data_;
}
inline bool Foam::pointEdgeStructuredWalk::valid() const
{
return previousPoint_ != vector::max;
}
// Checks for cyclic points
inline bool Foam::pointEdgeStructuredWalk::sameGeometry
(
const pointEdgeStructuredWalk& w2,
const scalar tol
) const
{
scalar diff = Foam::mag(dist() - w2.dist());
if (diff < SMALL)
{
return true;
}
else
{
if ((dist() > SMALL) && ((diff/dist()) < tol))
{
return true;
}
else
{
return false;
}
}
}
inline void Foam::pointEdgeStructuredWalk::leaveDomain
(
const polyPatch& patch,
const label patchPointI,
const point& coord
)
{
previousPoint_ -= coord;
}
inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
{
previousPoint_ = Foam::transform(rotTensor, previousPoint_);
}
// Update absolute geometric quantities. Note that distance (dist_)
// is not affected by leaving/entering domain.
inline void Foam::pointEdgeStructuredWalk::enterDomain
(
const polyPatch& patch,
const label patchPointI,
const point& coord
)
{
// back to absolute form
previousPoint_ += coord;
}
// Update this with information from connected edge
inline bool Foam::pointEdgeStructuredWalk::updatePoint
(
const polyMesh& mesh,
const label pointI,
const label edgeI,
const pointEdgeStructuredWalk& edgeInfo,
const scalar tol
)
{
if (inZone_)
{
return update(mesh.points()[pointI], edgeInfo, tol);
}
else
{
return false;
}
}
// Update this with new information on same point
inline bool Foam::pointEdgeStructuredWalk::updatePoint
(
const polyMesh& mesh,
const label pointI,
const pointEdgeStructuredWalk& newPointInfo,
const scalar tol
)
{
if (inZone_)
{
return update(mesh.points()[pointI], newPointInfo, tol);
}
else
{
return false;
}
}
// Update this with new information on same point. No extra information.
inline bool Foam::pointEdgeStructuredWalk::updatePoint
(
const pointEdgeStructuredWalk& newPointInfo,
const scalar tol
)
{
return update(newPointInfo, tol);
}
// Update this with information from connected point
inline bool Foam::pointEdgeStructuredWalk::updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointI,
const pointEdgeStructuredWalk& pointInfo,
const scalar tol
)
{
if (inZone_)
{
const pointField& points = mesh.points();
const edge& e = mesh.edges()[edgeI];
const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
return update(edgeMid, pointInfo, tol);
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline bool Foam::pointEdgeStructuredWalk::operator==
(
const Foam::pointEdgeStructuredWalk& rhs
) const
{
return previousPoint() == rhs.previousPoint();
}
inline bool Foam::pointEdgeStructuredWalk::operator!=
(
const Foam::pointEdgeStructuredWalk& rhs
) const
{
return !(*this == rhs);
}
// ************************************************************************* //

11
src/sampling/sampledSurface/distanceSurface/distanceSurface.C
View File

@ -30,8 +30,6 @@ License
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "isoSurface.H"
// #include "isoSurfaceCell.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -244,6 +242,14 @@ void Foam::distanceSurface::createGeometry()
distance_,
regularise_
)
//new isoSurfaceCell
//(
// fvm,
// cellDistance,
// pointDistance_,
// distance_,
// regularise_
//)
);
if (debug)
@ -284,6 +290,7 @@ Foam::distanceSurface::distanceSurface
distance_(readScalar(dict.lookup("distance"))),
signed_(readBool(dict.lookup("signed"))),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)),
zoneName_(word::null),
needsUpdate_(true),
isoSurfPtr_(NULL),

5
src/sampling/sampledSurface/distanceSurface/distanceSurface.H
View File

@ -68,6 +68,9 @@ class distanceSurface
//- Whether to coarsen
const Switch regularise_;
//- Whether to recalculate cell values as average of point values
const Switch average_;
//- zone name (if restricted to zones)
word zoneName_;
@ -82,6 +85,7 @@ class distanceSurface
scalarField pointDistance_;
//- Constructed iso surface
//autoPtr<isoSurfaceCell> isoSurfPtr_;
autoPtr<isoSurface> isoSurfPtr_;
//- triangles converted to faceList
@ -166,6 +170,7 @@ public:
}
//const isoSurfaceCell& surface() const
const isoSurface& surface() const
{
return isoSurfPtr_();

11
src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
View File

@ -25,7 +25,6 @@ License
\*---------------------------------------------------------------------------*/
#include "distanceSurface.H"
#include "isoSurface.H"
#include "volFieldsFwd.H"
#include "pointFields.H"
#include "volPointInterpolation.H"
@ -62,7 +61,15 @@ Foam::distanceSurface::interpolateField
);
// Sample.
return surface().interpolate(volFld, pointFld());
return surface().interpolate
(
(
average_
? pointAverage(pointFld())()
: volFld
),
pointFld()
);
}

103
src/sampling/sampledSurface/isoSurface/isoSurface.C
View File

@ -2054,61 +2054,68 @@ Foam::isoSurface::isoSurface
}
//if (false)
//{
List<FixedList<label, 3> > faceEdges;
labelList edgeFace0, edgeFace1;
Map<labelList> edgeFacesRest;
while (true)
if (false)
{
// Calculate addressing
calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
List<FixedList<label, 3> > faceEdges;
labelList edgeFace0, edgeFace1;
Map<labelList> edgeFacesRest;
// See if any dangling triangles
boolList keepTriangles;
label nDangling = markDanglingTriangles
(
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest,
keepTriangles
);
if (debug)
while (true)
{
Pout<< "isoSurface : detected " << nDangling
<< " dangling triangles." << endl;
}
if (nDangling == 0)
{
break;
}
// Create face map (new to old)
labelList subsetTriMap(findIndices(keepTriangles, true));
labelList subsetPointMap;
labelList reversePointMap;
triSurface::operator=
(
subsetMesh
// Calculate addressing
calcAddressing
(
*this,
subsetTriMap,
reversePointMap,
subsetPointMap
)
);
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
}
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest
);
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
//}
// See if any dangling triangles
boolList keepTriangles;
label nDangling = markDanglingTriangles
(
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest,
keepTriangles
);
if (debug)
{
Pout<< "isoSurface : detected " << nDangling
<< " dangling triangles." << endl;
}
if (nDangling == 0)
{
break;
}
// Create face map (new to old)
labelList subsetTriMap(findIndices(keepTriangles, true));
labelList subsetPointMap;
labelList reversePointMap;
triSurface::operator=
(
subsetMesh
(
*this,
subsetTriMap,
reversePointMap,
subsetPointMap
)
);
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
}
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
}
if (debug)

103
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
View File

@ -1398,6 +1398,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
)
:
mesh_(mesh),
cVals_(cVals),
pVals_(pVals),
iso_(iso),
mergeDistance_(mergeTol*mesh.bounds().mag())
{
@ -1563,59 +1565,68 @@ Foam::isoSurfaceCell::isoSurfaceCell
}
List<FixedList<label, 3> > faceEdges;
labelList edgeFace0, edgeFace1;
Map<labelList> edgeFacesRest;
while (true)
if (false)
{
// Calculate addressing
calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
List<FixedList<label, 3> > faceEdges;
labelList edgeFace0, edgeFace1;
Map<labelList> edgeFacesRest;
// See if any dangling triangles
boolList keepTriangles;
label nDangling = markDanglingTriangles
(
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest,
keepTriangles
);
if (debug)
while (true)
{
Pout<< "isoSurfaceCell : detected " << nDangling
<< " dangling triangles." << endl;
}
if (nDangling == 0)
{
break;
}
// Create face map (new to old)
labelList subsetTriMap(findIndices(keepTriangles, true));
labelList subsetPointMap;
labelList reversePointMap;
triSurface::operator=
(
subsetMesh
// Calculate addressing
calcAddressing
(
*this,
subsetTriMap,
reversePointMap,
subsetPointMap
)
);
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest
);
// See if any dangling triangles
boolList keepTriangles;
label nDangling = markDanglingTriangles
(
faceEdges,
edgeFace0,
edgeFace1,
edgeFacesRest,
keepTriangles
);
if (debug)
{
Pout<< "isoSurfaceCell : detected " << nDangling
<< " dangling triangles." << endl;
}
if (nDangling == 0)
{
break;
}
// Create face map (new to old)
labelList subsetTriMap(findIndices(keepTriangles, true));
labelList subsetPointMap;
labelList reversePointMap;
triSurface::operator=
(
subsetMesh
(
*this,
subsetTriMap,
reversePointMap,
subsetPointMap
)
);
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
}
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
}
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
//combineCellTriangles();
}

12
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
View File

@ -85,6 +85,10 @@ class isoSurfaceCell
//- Reference to mesh
const polyMesh& mesh_;
const scalarField& cVals_;
const scalarField& pVals_;
//- isoSurfaceCell value
const scalar iso_;
@ -352,6 +356,14 @@ public:
const Field<Type>& cCoords,
const Field<Type>& pCoords
) const;
//- Interpolates cCoords,pCoords.
template <class Type>
tmp<Field<Type> > interpolate
(
const Field<Type>& cCoords,
const Field<Type>& pCoords
) const;
};

52
src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
View File

@ -381,4 +381,56 @@ Foam::isoSurfaceCell::interpolate
}
template <class Type>
Foam::tmp<Foam::Field<Type> >
Foam::isoSurfaceCell::interpolate
(
const Field<Type>& cCoords,
const Field<Type>& pCoords
) const
{
DynamicList<Type> triPoints(nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data
DynamicList<Type> snappedPoints;
labelList snappedCc(mesh_.nCells(), -1);
labelList snappedPoint(mesh_.nPoints(), -1);
generateTriPoints
(
cVals_,
pVals_,
cCoords,
pCoords,
snappedPoints,
snappedCc,
snappedPoint,
triPoints,
triMeshCells
);
// One value per point
tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
Field<Type>& values = tvalues();
forAll(triPoints, i)
{
label mergedPointI = triPointMergeMap_[i];
if (mergedPointI >= 0)
{
values[mergedPointI] = triPoints[i];
}
}
return tvalues;
}
// ************************************************************************* //

100
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
View File

@ -164,7 +164,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
// point field.
if (average_)
{
storedVolFieldPtr_.reset(average(fvm, *pointFieldPtr_).ptr());
storedVolFieldPtr_.reset
(
pointAverage(*pointFieldPtr_).ptr()
);
volFieldPtr_ = storedVolFieldPtr_.operator->();
}
@ -265,7 +268,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
{
storedVolSubFieldPtr_.reset
(
average(subFvm, *pointSubFieldPtr_).ptr()
pointAverage(*pointSubFieldPtr_).ptr()
);
volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
}
@ -286,99 +289,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
}
Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
(
const fvMesh& mesh,
const pointScalarField& pfld
) const
{
tmp<volScalarField> tcellAvg
(
new volScalarField
(
IOobject
(
"cellAvg",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, scalar(0.0))
)
);
volScalarField& cellAvg = tcellAvg();
labelField nPointCells(mesh.nCells(), 0);
{
for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
{
const labelList& pCells = mesh.pointCells(pointI);
forAll(pCells, i)
{
label cellI = pCells[i];
cellAvg[cellI] += pfld[pointI];
nPointCells[cellI]++;
}
}
}
forAll(cellAvg, cellI)
{
cellAvg[cellI] /= nPointCells[cellI];
}
// Give value to calculatedFvPatchFields
cellAvg.correctBoundaryConditions();
return tcellAvg;
}
Foam::tmp<Foam::pointScalarField> Foam::sampledIsoSurface::average
(
const pointMesh& pMesh,
const volScalarField& fld
) const
{
tmp<pointScalarField> tpointAvg
(
new pointScalarField
(
IOobject
(
"pointAvg",
fld.time().timeName(),
fld.db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
pMesh,
dimensionedScalar("zero", dimless, scalar(0.0))
)
);
pointScalarField& pointAvg = tpointAvg();
for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
{
const labelList& pCells = fld.mesh().pointCells(pointI);
forAll(pCells, i)
{
pointAvg[pointI] += fld[pCells[i]];
}
pointAvg[pointI] /= pCells.size();
}
// Give value to calculatedFvPatchFields
pointAvg.correctBoundaryConditions();
return tpointAvg;
}
bool Foam::sampledIsoSurface::updateGeometry() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());

12
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
View File

@ -118,18 +118,6 @@ class sampledIsoSurface
//- Get fields needed to recreate iso surface.
void getIsoFields() const;
tmp<volScalarField> average
(
const fvMesh&,
const pointScalarField&
) const;
tmp<pointScalarField> average
(
const pointMesh&,
const volScalarField& fld
) const;
//- Create iso surface (if time has changed)
// Do nothing (and return false) if no update was needed
bool updateGeometry() const;

20
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
View File

@ -71,7 +71,15 @@ Foam::sampledIsoSurface::interpolateField
volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
// Sample.
return surface().interpolate(volSubFld, tpointSubFld());
return surface().interpolate
(
(
average_
? pointAverage(tpointSubFld())()
: volSubFld
),
tpointSubFld()
);
}
else
{
@ -79,7 +87,15 @@ Foam::sampledIsoSurface::interpolateField
volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
// Sample.
return surface().interpolate(volFld, tpointFld());
return surface().interpolate
(
(
average_
? pointAverage(tpointFld())()
: volFld
),
tpointFld()
);
}
}

12
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
View File

@ -30,7 +30,6 @@ License
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "isoSurface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -173,7 +172,7 @@ void Foam::sampledCuttingPlane::createGeometry()
forAll(fld, i)
{
fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
}
}
}
@ -231,6 +230,14 @@ void Foam::sampledCuttingPlane::createGeometry()
0.0,
regularise_
)
//new isoSurfaceCell
//(
// fvm,
// cellDistance,
// pointDistance_,
// 0.0,
// regularise_
//)
);
if (debug)
@ -254,6 +261,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
plane_(dict),
mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)),
zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
exposedPatchName_(word::null),
needsUpdate_(true),

6
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
View File

@ -38,6 +38,7 @@ SourceFiles
#include "sampledSurface.H"
#include "isoSurface.H"
//#include "isoSurfaceCell.H"
#include "plane.H"
#include "ZoneIDs.H"
#include "fvMeshSubset.H"
@ -66,6 +67,9 @@ class sampledCuttingPlane
//- Whether to coarsen
const Switch regularise_;
//- Whether to recalculate cell values as average of point values
const Switch average_;
//- zone name/index (if restricted to zones)
mutable cellZoneID zoneID_;
@ -86,6 +90,7 @@ class sampledCuttingPlane
scalarField pointDistance_;
//- Constructed iso surface
//autoPtr<isoSurfaceCell> isoSurfPtr_;
autoPtr<isoSurface> isoSurfPtr_;
//- triangles converted to faceList
@ -170,6 +175,7 @@ public:
}
//const isoSurfaceCell& surface() const
const isoSurface& surface() const
{
return isoSurfPtr_();

21
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
View File

@ -26,7 +26,6 @@ License
#include "volPointInterpolation.H"
#include "sampledCuttingPlane.H"
#include "isoSurface.H"
#include "volFieldsFwd.H"
#include "pointFields.H"
#include "volPointInterpolation.H"
@ -67,7 +66,15 @@ Foam::sampledCuttingPlane::interpolateField
volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
// Sample.
return surface().interpolate(volSubFld, tpointSubFld());
return surface().interpolate
(
(
average_
? pointAverage(tpointSubFld())()
: volSubFld
),
tpointSubFld()
);
}
else
{
@ -77,7 +84,15 @@ Foam::sampledCuttingPlane::interpolateField
);
// Sample.
return surface().interpolate(volFld, tpointFld());
return surface().interpolate
(
(
average_
? pointAverage(tpointFld())()
: volFld
),
tpointFld()
);
}
}

7
src/sampling/sampledSurface/sampledSurface/sampledSurface.H
View File

@ -312,6 +312,13 @@ public:
//- Project field onto surface
tmp<Field<vector> > project(const Field<tensor>&) const;
//- Interpolate from points to cell centre
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > pointAverage
(
const GeometricField<Type, pointPatchField, pointMesh>& pfld
) const;
//- Sample field on surface
virtual tmp<scalarField> sample
(

55
src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
View File

@ -155,4 +155,59 @@ Foam::sampledSurface::project
}
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::sampledSurface::pointAverage
(
const GeometricField<Type, pointPatchField, pointMesh>& pfld
) const
{
const fvMesh& mesh = dynamic_cast<const fvMesh&>(pfld.mesh()());
tmp<GeometricField<Type, fvPatchField, volMesh> > tcellAvg
(
new GeometricField<Type, fvPatchField, volMesh>
(
IOobject
(
"cellAvg",
mesh.time().timeName(),
pfld.db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
)
);
GeometricField<Type, fvPatchField, volMesh>& cellAvg = tcellAvg();
labelField nPointCells(mesh.nCells(), 0);
{
for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
{
const labelList& pCells = mesh.pointCells(pointI);
forAll(pCells, i)
{
label cellI = pCells[i];
cellAvg[cellI] += pfld[pointI];
nPointCells[cellI]++;
}
}
}
forAll(cellAvg, cellI)
{
cellAvg[cellI] /= nPointCells[cellI];
}
// Give value to calculatedFvPatchFields
cellAvg.correctBoundaryConditions();
return tcellAvg;
}
// ************************************************************************* //

209
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,11 +25,10 @@ License
\*---------------------------------------------------------------------------*/
#include "sampledTriSurfaceMesh.H"
#include "dictionary.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "volFields.H"
#include "treeDataPoint.H"
#include "meshSearch.H"
#include "Tuple2.H"
#include "globalIndex.H"
#include "addToRunTimeSelectionTable.H"
@ -44,6 +43,25 @@ namespace Foam
sampledTriSurfaceMesh,
word
);
//- Private class for finding nearest
// - global index
// - sqr(distance)
typedef Tuple2<scalar, label> nearInfo;
class nearestEqOp
{
public:
void operator()(nearInfo& x, const nearInfo& y) const
{
if (y.first() < x.first())
{
x = y;
}
}
};
}
@ -63,16 +81,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
(
name,
mesh.time().constant(), // instance
"triSurface", // local
mesh, // registry
"triSurface", // local
mesh, // registry
IOobject::MUST_READ,
IOobject::NO_WRITE
IOobject::NO_WRITE,
false
)
),
needsUpdate_(true),
cellLabels_(0),
pointMap_(0),
faceMap_(0)
pointToFace_(0)
{}
@ -90,16 +108,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
(
dict.lookup("surface"),
mesh.time().constant(), // instance
"triSurface", // local
mesh, // registry
"triSurface", // local
mesh, // registry
IOobject::MUST_READ,
IOobject::NO_WRITE
IOobject::NO_WRITE,
false
)
),
needsUpdate_(true),
cellLabels_(0),
pointMap_(0),
faceMap_(0)
pointToFace_(0)
{}
@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
sampledSurface::clearGeom();
MeshStorage::clear();
cellLabels_.clear();
pointMap_.clear();
faceMap_.clear();
pointToFace_.clear();
needsUpdate_ = true;
return true;
@ -146,95 +163,139 @@ bool Foam::sampledTriSurfaceMesh::update()
// Find the cells the triangles of the surface are in.
// Does approximation by looking at the face centres only
const pointField fc = surface_.faceCentres();
cellLabels_.setSize(fc.size());
const pointField& fc = surface_.faceCentres();
meshSearch meshSearcher(mesh(), false);
const indexedOctree<treeDataPoint>& cellCentreTree =
meshSearcher.cellCentreTree();
// Global numbering for cells - only used to uniquely identify local cells.
globalIndex globalCells(mesh().nCells());
List<nearInfo> nearest(fc.size());
forAll(nearest, i)
{
nearest[i].first() = GREAT;
nearest[i].second() = labelMax;
}
// Search triangles using nearest on local mesh
forAll(fc, triI)
{
cellLabels_[triI] = meshSearcher.findCell
pointIndexHit nearInfo = cellCentreTree.findNearest
(
fc[triI],
-1, // seed
true // use tree
sqr(GREAT)
);
}
boolList include(surface_.size(), true);
label nNotFound = 0;
forAll(cellLabels_, triI)
{
if (cellLabels_[triI] == -1)
if (nearInfo.hit())
{
include[triI] = false;
nNotFound++;
nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
}
}
label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
// See which processor has the nearest.
Pstream::listCombineGather(nearest, nearestEqOp());
Pstream::listCombineScatter(nearest);
boolList include(surface_.size(), false);
cellLabels_.setSize(fc.size());
cellLabels_ = -1;
label nFound = 0;
forAll(nearest, triI)
{
if (nearest[triI].second() == labelMax)
{
// Not found on any processor. How to map?
}
else if (globalCells.isLocal(nearest[triI].second()))
{
cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
include[triI] = true;
nFound++;
}
}
if (debug)
{
Pout<< "surface:" << surface_.size()
<< " included:" << surface_.size()-nNotFound
<< " total not found" << nTotalNotFound << endl;
Pout<< "Local out of faces:" << cellLabels_.size()
<< " keeping:" << nFound << endl;
}
// Now subset the surface. Do not use triSurface::subsetMesh since requires
// original surface to be in compact numbering.
const triSurface& s = surface_;
// Compact to original triangle
labelList faceMap(s.size());
// Compact to original points
labelList pointMap(s.points().size());
// From original point to compact points
labelList reversePointMap(s.points().size(), -1);
// Add to master all triangles are outside all meshes.
{
boolList onAnyProc(surface_.size(), false);
Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
label newPointI = 0;
label newTriI = 0;
if (Pstream::master())
forAll(s, triI)
{
label nAdded = 0;
forAll(onAnyProc, triI)
if (include[triI])
{
if (!onAnyProc[triI])
faceMap[newTriI++] = triI;
const labelledTri& f = s[triI];
forAll(f, fp)
{
include[triI] = true;
nAdded++;
if (reversePointMap[f[fp]] == -1)
{
pointMap[newPointI] = f[fp];
reversePointMap[f[fp]] = newPointI++;
}
}
}
}
faceMap.setSize(newTriI);
pointMap.setSize(newPointI);
}
if (debug)
{
Pout<< "nAdded to master:" << nAdded << endl;
}
// Subset cellLabels
cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
// Store any face per point
pointToFace_.setSize(pointMap.size());
// Create faces and points for subsetted surface
faceList& faces = this->storedFaces();
faces.setSize(faceMap.size());
forAll(faceMap, i)
{
const triFace& f = s[faceMap[i]];
triFace newF
(
reversePointMap[f[0]],
reversePointMap[f[1]],
reversePointMap[f[2]]
);
faces[i] = newF.triFaceFace();
forAll(newF, fp)
{
pointToFace_[newF[fp]] = i;
}
}
// Now subset the surface
triSurface localSurface
(
surface_.subsetMesh
(
include,
pointMap_,
faceMap_
)
);
// And convert into faces.
faceList& faces = this->storedFaces();
faces.setSize(localSurface.size());
forAll(localSurface, i)
{
faces[i] = localSurface[i].triFaceFace();
}
this->storedPoints() = localSurface.points();
this->storedPoints() = pointField(s.points(), pointMap);
if (debug)
{
print(Pout);
Pout << endl;
Pout<< endl;
}
needsUpdate_ = false;

17
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,15 +29,13 @@ Description
A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
of the triSurface.
It samples using the cell nearest to the triangle centre so does
not check the cell the centre is actually in ...
In parallel every processor just operates on the part of the surface
where the face centres are inside the mesh. It is then up to the
caller to stitch the partial surfaces together.
No check is done for triangle centres being on multiple processors
(should never happen but truncation errors ...)
Any value on a triangle outside the mesh will be set to 0.
SourceFiles
sampledTriSurfaceMesh.C
@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
//- Track if the surface needs an update
mutable bool needsUpdate_;
//- Local cell labels
//- From local surface triangle to mesh cell.
labelList cellLabels_;
//- From local surface back to surface_
labelList pointMap_;
//- From local surface back to surface_
labelList faceMap_;
labelList pointToFace_;
// Private Member Functions

33
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
View File

@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
) const
{
// One value per face
tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
Field<Type>& values = tvalues();
forAll(faceMap_, i)
forAll(cellLabels_, triI)
{
label cellI = cellLabels_[faceMap_[i]];
if (cellI != -1)
{
values[i] = vField[cellI];
}
else
{
values[i] = pTraits<Type>::zero;
}
values[triI] = vField[cellLabels_[triI]];
}
return tvalues;
@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
) const
{
// One value per vertex
tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
Field<Type>& values = tvalues();
forAll(pointMap_, i)
forAll(pointToFace_, pointI)
{
label origPointI = pointMap_[i];
label origTriI = surface_.pointFaces()[origPointI][0];
label triI = pointToFace_[pointI];
label cellI = cellLabels_[triI];
label cellI = cellLabels_[origTriI];
if (cellI != -1)
{
values[i] = interpolator.interpolate(points()[i], cellI);
}
else
{
values[i] = pTraits<Type>::zero;
}
values[pointI] = interpolator.interpolate(points()[pointI], cellI);
}
return tvalues;