Merge branch 'master' of ssh://opencfd:8007/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2013-04-27 15:31:26 +01:00
74 changed files with 5378 additions and 1932 deletions

View File

@ -17,7 +17,7 @@ rAU = 1.0/UEqn().A();
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p) + fvOptions(U));
solve(UEqn() == -fvc::grad(p));
fvOptions.correct(U);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,6 +47,8 @@ Description
#include "snapParameters.H"
#include "layerParameters.H"
#include "vtkSetWriter.H"
#include "faceSet.H"
#include "motionSmoother.H"
using namespace Foam;
@ -675,6 +677,31 @@ int main(int argc, char *argv[])
}
{
// Check final mesh
Info<< "Checking final mesh ..." << endl;
faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
const label nErrors = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
if (nErrors > 0)
{
Info<< "Finished meshing with " << nErrors << " illegal faces"
<< " (concave, zero area or negative cell pyramid volume)"
<< endl;
wrongFaces.write();
}
else
{
Info<< "Finished meshing without any errors" << endl;
}
}
Info<< "Finished meshing in = "
<< runTime.elapsedCpuTime() << " s." << endl;

View File

@ -832,5 +832,39 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
if (allGeometry)
{
faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
if (mesh.checkFaceWeight(true, 0.05, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with low interpolation weights to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
if (allGeometry)
{
faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
if (mesh.checkVolRatio(true, 0.05, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with low volume ratio cells to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
return noFailedChecks;
}

View File

@ -138,6 +138,8 @@ baffles
//- Select faces and orientation through a searchableSurface
type searchableSurface;
surface searchablePlate;
//name sphere.stl; // name if surface=triSurfaceMesh
origin (0.099 -0.006 0.004);
span (0 0.012 0.012);

View File

@ -28,6 +28,7 @@ License
#include "syncTools.H"
#include "searchableSurface.H"
#include "fvMesh.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -61,7 +62,15 @@ Foam::faceSelections::searchableSurfaceSelection::searchableSurfaceSelection
searchableSurface::New
(
word(dict.lookup("surface")),
IOobject
(
dict.lookupOrDefault("name", mesh.objectRegistry::db().name()),
mesh.time().constant(),
"triSurface",
mesh.objectRegistry::db(),
IOobject::MUST_READ,
IOobject::NO_WRITE
),
dict
)
)

View File

@ -378,6 +378,7 @@ FoamFile
// surface searchableSphere;
// centre (0.05 0.05 0.005);
// radius 0.025;
// //name sphere.stl; // Optional name if surface triSurfaceMesh
// }
// }
//

View File

@ -2,8 +2,10 @@
cd ${0%/*} || exit 1 # run from this directory
#set -x
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
if [ "$ParaView_VERSION" == "3.98.1" ]
then
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
[ -n "$PV_PLUGIN_PATH" ] || {
echo "$0 : PV_PLUGIN_PATH not valid - it is unset"
exit 1
@ -15,8 +17,11 @@ then
wmake libso vtkPV398Readers
PV398blockMeshReader/Allwmake
PV398FoamReader/Allwmake
else
else
echo "ERROR: ParaView not found in $ParaView_DIR"
fi
else
echo "WARN: PV398 readers not building: ParaView_VERSION=$ParaView_VERSION"
fi
# ----------------------------------------------------------------- end-of-file

View File

@ -2,8 +2,10 @@
cd ${0%/*} || exit 1 # run from this directory
#set -x
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
if [ "$ParaView_VERSION" != "3.98.1" ]
then
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
[ -n "$PV_PLUGIN_PATH" ] || {
echo "$0 : PV_PLUGIN_PATH not valid - it is unset"
exit 1
@ -15,8 +17,11 @@ then
wmake libso vtkPV3Readers
PV3blockMeshReader/Allwmake
PV3FoamReader/Allwmake
else
else
echo "ERROR: ParaView not found in $ParaView_DIR"
fi
else
echo "WARN: PV3 readers not building: ParaView_VERSION=$ParaView_VERSION"
fi
# ----------------------------------------------------------------- end-of-file

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -179,7 +179,10 @@ void Foam::vtkPV3Foam::convertMeshPatches
const word patchName = getPartName(partId);
labelHashSet patchIds(patches.patchSet(List<wordRe>(1, patchName)));
labelHashSet patchIds
(
patches.patchSet(List<wordRe>(1, wordRe(patchName)))
);
if (debug)
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -182,8 +182,17 @@ void Foam::csvTableReader<Type>::write(Ostream& os) const
<< headerLine_ << token::END_STATEMENT << nl;
os.writeKeyword("timeColumn")
<< timeColumn_ << token::END_STATEMENT << nl;
os.writeKeyword("valueColumns")
<< componentColumns_ << token::END_STATEMENT << nl;
// Force writing labelList in ascii
os.writeKeyword("valueColumns");
if (os.format() == IOstream::BINARY)
{
os.format(IOstream::ASCII);
os << componentColumns_;
os.format(IOstream::BINARY);
}
os << token::END_STATEMENT << nl;
os.writeKeyword("separator")
<< string(separator_) << token::END_STATEMENT << nl;
}

View File

@ -265,6 +265,24 @@ private:
const Vector<label>& meshD
) const;
bool checkFaceWeight
(
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const scalar minWeight,
labelHashSet* setPtr
) const;
bool checkVolRatio
(
const scalarField& cellVols,
const bool report,
const scalar minRatio,
labelHashSet* setPtr
) const;
public:
// Public typedefs
@ -612,6 +630,22 @@ public:
const bool detailedReport = false
) const;
//- Check for face weights
virtual bool checkFaceWeight
(
const bool report,
const scalar minWeight = 0.05,
labelHashSet* setPtr = NULL
) const;
//- Check for neighbouring cell volumes
virtual bool checkVolRatio
(
const bool report,
const scalar minRatio = 0.01,
labelHashSet* setPtr = NULL
) const;
// Helper functions

View File

@ -142,7 +142,8 @@ bool Foam::polyMesh::checkFaceOrthogonality
if (severeNonOrth > 0)
{
Info<< " *Number of severely non-orthogonal faces: "
Info<< " *Number of severely non-orthogonal (> "
<< primitiveMesh::nonOrthThreshold_ << " degrees) faces: "
<< severeNonOrth << "." << endl;
}
}
@ -427,6 +428,8 @@ bool Foam::polyMesh::checkCellDeterminant
const Vector<label>& meshD
) const
{
const scalar warnDet = 1e-3;
if (debug)
{
Info<< "bool polyMesh::checkCellDeterminant(const bool"
@ -450,7 +453,7 @@ bool Foam::polyMesh::checkCellDeterminant
forAll (cellDeterminant, cellI)
{
if (cellDeterminant[cellI] < 1e-3)
if (cellDeterminant[cellI] < warnDet)
{
if (setPtr)
{
@ -480,7 +483,8 @@ bool Foam::polyMesh::checkCellDeterminant
{
if (debug || report)
{
Info<< " ***Cells with small determinant found, number of cells: "
Info<< " ***Cells with small determinant (< "
<< warnDet << ") found, number of cells: "
<< nErrorCells << endl;
}
@ -500,6 +504,192 @@ bool Foam::polyMesh::checkCellDeterminant
}
bool Foam::polyMesh::checkFaceWeight
(
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const scalar minWeight,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool polyMesh::checkFaceWeight(const bool"
<< ", labelHashSet*) const: "
<< "checking for low face interpolation weights" << endl;
}
tmp<scalarField> tfaceWght = polyMeshTools::faceWeights
(
*this,
fCtrs,
fAreas,
cellCtrs
);
scalarField& faceWght = tfaceWght();
label nErrorFaces = 0;
scalar minDet = GREAT;
scalar sumDet = 0.0;
label nSummed = 0;
// Statistics only for internal and masters of coupled faces
PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
forAll (faceWght, faceI)
{
if (faceWght[faceI] < minWeight)
{
// Note: insert both sides of coupled faces
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorFaces++;
}
// Note: statistics only on master of coupled faces
if (isMasterFace[faceI])
{
minDet = min(minDet, faceWght[faceI]);
sumDet += faceWght[faceI];
nSummed++;
}
}
reduce(nErrorFaces, sumOp<label>());
reduce(minDet, minOp<scalar>());
reduce(sumDet, sumOp<scalar>());
reduce(nSummed, sumOp<label>());
if (debug || report)
{
if (nSummed > 0)
{
Info<< " Face interpolation weight : minimum: " << minDet
<< " average: " << sumDet/nSummed
<< endl;
}
}
if (nErrorFaces > 0)
{
if (debug || report)
{
Info<< " ***Faces with small interpolation weight (< " << minWeight
<< ") found, number of faces: "
<< nErrorFaces << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Face interpolation weight check OK." << endl;
}
return false;
}
return false;
}
bool Foam::polyMesh::checkVolRatio
(
const scalarField& cellVols,
const bool report,
const scalar minRatio,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool polyMesh::checkVolRatio(const bool"
<< ", labelHashSet*) const: "
<< "checking for volume ratio < " << minRatio << endl;
}
tmp<scalarField> tvolRatio = polyMeshTools::volRatio(*this, cellVols);
scalarField& volRatio = tvolRatio();
label nErrorFaces = 0;
scalar minDet = GREAT;
scalar sumDet = 0.0;
label nSummed = 0;
// Statistics only for internal and masters of coupled faces
PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
forAll (volRatio, faceI)
{
if (volRatio[faceI] < minRatio)
{
// Note: insert both sides of coupled faces
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorFaces++;
}
// Note: statistics only on master of coupled faces
if (isMasterFace[faceI])
{
minDet = min(minDet, volRatio[faceI]);
sumDet += volRatio[faceI];
nSummed++;
}
}
reduce(nErrorFaces, sumOp<label>());
reduce(minDet, minOp<scalar>());
reduce(sumDet, sumOp<scalar>());
reduce(nSummed, sumOp<label>());
if (debug || report)
{
if (nSummed > 0)
{
Info<< " Face volume ratio : minimum: " << minDet
<< " average: " << sumDet/nSummed
<< endl;
}
}
if (nErrorFaces > 0)
{
if (debug || report)
{
Info<< " ***Faces with small volume ratio (< " << minRatio
<< ") found, number of faces: "
<< nErrorFaces << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Face volume ratio check OK." << endl;
}
return false;
}
return false;
}
//- Could override checkClosedBoundary to not look at (collocated!) coupled
// faces
//bool Foam::polyMesh::checkClosedBoundary(const bool report) const
@ -582,6 +772,36 @@ bool Foam::polyMesh::checkCellDeterminant
}
bool Foam::polyMesh::checkFaceWeight
(
const bool report,
const scalar minWeight,
labelHashSet* setPtr
) const
{
return checkFaceWeight
(
faceCentres(),
faceAreas(),
cellCentres(),
report,
minWeight,
setPtr
);
}
bool Foam::polyMesh::checkVolRatio
(
const bool report,
const scalar minRatio,
labelHashSet* setPtr
) const
{
return checkVolRatio(cellVolumes(), report, minRatio, setPtr);
}
bool Foam::polyMesh::checkMeshMotion
(
const pointField& newPoints,

View File

@ -187,4 +187,112 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
}
Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
scalarField& weight = tweight();
// Internal faces
forAll(nei, faceI)
{
const point& fc = fCtrs[faceI];
const vector& fa = fAreas[faceI];
scalar dOwn = mag(fa & (fc-cellCtrs[own[faceI]]));
scalar dNei = mag(fa & (cellCtrs[nei[faceI]]-fc));
weight[faceI] = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
}
// Coupled faces
pointField neiCc;
syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start() + i;
label bFaceI = faceI - mesh.nInternalFaces();
const point& fc = fCtrs[faceI];
const vector& fa = fAreas[faceI];
scalar dOwn = mag(fa & (fc-cellCtrs[own[faceI]]));
scalar dNei = mag(fa & (neiCc[bFaceI]-fc));
weight[faceI] = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
}
}
}
return tweight;
}
Foam::tmp<Foam::scalarField> Foam::polyMeshTools::volRatio
(
const polyMesh& mesh,
const scalarField& vol
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
scalarField& ratio = tratio();
// Internal faces
forAll(nei, faceI)
{
scalar volOwn = vol[own[faceI]];
scalar volNei = vol[nei[faceI]];
ratio[faceI] = min(volOwn,volNei)/(max(volOwn, volNei)+VSMALL);
}
// Coupled faces
scalarField neiVol;
syncTools::swapBoundaryCellList(mesh, vol, neiVol);
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start() + i;
label bFaceI = faceI - mesh.nInternalFaces();
scalar volOwn = vol[own[faceI]];
scalar volNei = neiVol[bFaceI];
ratio[faceI] = min(volOwn,volNei)/(max(volOwn, volNei)+VSMALL);
}
}
}
return tratio;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -95,6 +95,22 @@ public:
// );
// }
//- Generate interpolation factors field
static tmp<scalarField> faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate volume ratio field
static tmp<scalarField> volRatio
(
const polyMesh& mesh,
const scalarField& vol
);
};

View File

@ -63,8 +63,6 @@ class processorCyclicPolyPatch
//- Index of originating patch
mutable label referPatchID_;
// Private member functions
protected:
@ -223,7 +221,6 @@ public:
// Destructor
virtual ~processorCyclicPolyPatch();

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -42,6 +42,7 @@ SourceFiles
#include "Random.H"
#include "FixedList.H"
#include "UList.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -222,6 +223,15 @@ public:
label& nearLabel
) const;
//- Return nearest point to line on triangle. Returns hit if
// point is inside triangle. Sets edgePoint to point on edge
// (hit if nearest is inside line)
inline pointHit nearestPoint
(
const linePointRef& edge,
pointHit& edgePoint
) const;
// IOstream operators

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -659,6 +659,132 @@ inline bool Foam::triangle<Point, PointRef>::classify
}
template<class Point, class PointRef>
inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
(
const linePointRef& ln,
pointHit& lnInfo
) const
{
vector q = ln.vec();
pointHit triInfo
(
triangle<Point, PointRef>::intersection
(
ln.start(),
q,
intersection::FULL_RAY
)
);
if (triInfo.hit())
{
// Line hits triangle. Find point on line.
if (triInfo.distance() > 1)
{
// Hit beyond endpoint
lnInfo.setMiss(true);
lnInfo.setPoint(ln.end());
scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
lnInfo.setDistance(dist);
triInfo.setMiss(true);
triInfo.setDistance(dist);
}
else if (triInfo.distance() < 0)
{
// Hit beyond startpoint
lnInfo.setMiss(true);
lnInfo.setPoint(ln.start());
scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
lnInfo.setDistance(dist);
triInfo.setMiss(true);
triInfo.setDistance(dist);
}
else
{
// Hit on line
lnInfo.setHit();
lnInfo.setPoint(triInfo.hitPoint());
lnInfo.setDistance(0.0);
triInfo.setDistance(0.0);
}
}
else
{
// Line skips triangle. See which triangle edge it gets closest to
point nearestEdgePoint;
point nearestLinePoint;
label minEdgeIndex = 0;
scalar minDist = ln.nearestDist
(
linePointRef(a_, b_),
nearestLinePoint,
nearestEdgePoint
);
{
point linePoint;
point triEdgePoint;
scalar dist = ln.nearestDist
(
linePointRef(b_, c_),
linePoint,
triEdgePoint
);
if (dist < minDist)
{
minDist = dist;
nearestEdgePoint = triEdgePoint;
nearestLinePoint = linePoint;
minEdgeIndex = 1;
}
}
{
point linePoint;
point triEdgePoint;
scalar dist = ln.nearestDist
(
linePointRef(c_, a_),
linePoint,
triEdgePoint
);
if (dist < minDist)
{
minDist = dist;
nearestEdgePoint = triEdgePoint;
nearestLinePoint = linePoint;
minEdgeIndex = 2;
}
}
lnInfo.setDistance(minDist);
triInfo.setDistance(minDist);
triInfo.setMiss(false);
triInfo.setPoint(nearestEdgePoint);
// Convert point on line to pointHit
if (Foam::mag(nearestLinePoint-ln.start()) < SMALL)
{
lnInfo.setMiss(true);
lnInfo.setPoint(ln.start());
}
else if (Foam::mag(nearestLinePoint-ln.end()) < SMALL)
{
lnInfo.setMiss(true);
lnInfo.setPoint(ln.end());
}
else
{
lnInfo.setHit();
lnInfo.setPoint(nearestLinePoint);
}
}
return triInfo;
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Point, class PointRef>

View File

@ -74,8 +74,17 @@ void Foam::CSV<Type>::writeData(Ostream& os) const
os.writeKeyword("nHeaderLine") << nHeaderLine_ << token::END_STATEMENT
<< nl;
os.writeKeyword("refColumn") << refColumn_ << token::END_STATEMENT << nl;
os.writeKeyword("componentColumns") << componentColumns_
<< token::END_STATEMENT << nl;
// Force writing labelList in ascii
os.writeKeyword("componentColumns");
if (os.format() == IOstream::BINARY)
{
os.format(IOstream::ASCII);
os << componentColumns_;
os.format(IOstream::BINARY);
}
os << token::END_STATEMENT << nl;
os.writeKeyword("separator") << string(separator_)
<< token::END_STATEMENT << nl;
os.writeKeyword("mergeSeparators") << string(mergeSeparators_)

View File

@ -120,22 +120,25 @@ public:
inline wordRe(const wordRe&);
//- Construct from keyType
inline wordRe(const keyType&, const compOption=LITERAL);
inline explicit wordRe(const keyType&);
//- Construct from keyType
inline wordRe(const keyType&, const compOption);
//- Construct as copy of word
inline explicit wordRe(const word&);
//- Construct as copy of character array
// Optionally specify how it should be treated.
inline wordRe(const char*, const compOption = LITERAL);
inline explicit wordRe(const char*, const compOption = LITERAL);
//- Construct as copy of string.
// Optionally specify how it should be treated.
inline wordRe(const string&, const compOption = LITERAL);
inline explicit wordRe(const string&, const compOption = LITERAL);
//- Construct as copy of std::string
// Optionally specify how it should be treated.
inline wordRe(const std::string&, const compOption = LITERAL);
inline explicit wordRe(const std::string&, const compOption = LITERAL);
//- Construct from Istream
// Words are treated as literals, strings with an auto-test

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -65,6 +65,18 @@ inline Foam::wordRe::wordRe(const word& str)
{}
inline Foam::wordRe::wordRe(const keyType& str)
:
word(str, false),
re_()
{
if (str.isPattern())
{
compile();
}
}
inline Foam::wordRe::wordRe(const keyType& str, const compOption opt)
:
word(str, false),

View File

@ -678,34 +678,8 @@ void Foam::motionSmoother::correct()
}
void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
void Foam::motionSmoother::setDisplacementPatchFields()
{
// See comment in .H file about shared points.
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
const labelList& meshPoints = pp.meshPoints();
forAll(meshPoints, i)
{
displacement_[meshPoints[i]] = vector::zero;
}
}
}
const labelList& ppMeshPoints = pp_.meshPoints();
// Set internal point data from displacement on combined patch points.
forAll(ppMeshPoints, patchPointI)
{
displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
}
// Adapt the fixedValue bc's (i.e. copy internal point data to
// boundaryField for all affected patches)
forAll(adaptPatchIDs_, i)
@ -765,6 +739,42 @@ void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
displacement_.boundaryField()[patchI] ==
displacement_.boundaryField()[patchI].patchInternalField();
}
}
void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
{
// See comment in .H file about shared points.
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
const labelList& meshPoints = pp.meshPoints();
forAll(meshPoints, i)
{
displacement_[meshPoints[i]] = vector::zero;
}
}
}
const labelList& ppMeshPoints = pp_.meshPoints();
// Set internal point data from displacement on combined patch points.
forAll(ppMeshPoints, patchPointI)
{
displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
}
// Adapt the fixedValue bc's (i.e. copy internal point data to
// boundaryField for all affected patches)
setDisplacementPatchFields();
if (debug)
{

View File

@ -385,6 +385,10 @@ public:
//- Take over existing mesh position.
void correct();
//- Set patch fields on displacement to be consistent with
// internal values.
void setDisplacementPatchFields();
//- Set displacement field from displacement on patch points.
// Modify provided displacement to be consistent with actual
// boundary conditions on displacement. Note: resets the

View File

@ -25,7 +25,18 @@ Class
Foam::mappedPatchFieldBase
Description
Functionality for sampling fields using mappedPatchBase.
Functionality for sampling fields using mappedPatchBase. Every call to
mappedField() returns a sampled field, optionally scaled to maintain an
area-weighted average.
Example usage:
{
fieldName T; // default is same as fvPatchField
setAverage false;
average 1.0; // only if setAverage=true
interpolationScheme cellPoint; // default is cell
}
SourceFiles
mappedPatchFieldBase.C

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,7 +29,7 @@ Group
Description
This velocity inlet/outlet boundary condition is applied to pressure
boundaries where the pressure is specified. A zero-gradient condtion is
boundaries where the pressure is specified. A zero-gradient condition is
applied for outflow (as defined by the flux); for inflow, the velocity is
obtained from the patch-face normal component of the internal-cell value.

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -62,7 +62,11 @@ Foam::fv::fourthGrad<Type>::calcGrad
// Assemble the second-order least-square gradient
// Calculate the second-order least-square gradient
tmp<GeometricField<GradType, fvPatchField, volMesh> > tsecondfGrad
= leastSquaresGrad<Type>(mesh).grad(vsf);
= leastSquaresGrad<Type>(mesh).grad
(
vsf,
"leastSquaresGrad(" + vsf.name() + ")"
);
const GeometricField<GradType, fvPatchField, volMesh>& secondfGrad =
tsecondfGrad();

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -477,7 +477,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
}
}
}

View File

@ -466,7 +466,7 @@ void Foam::ReactingParcel<ParcelType>::calc
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -260,7 +260,7 @@ void Foam::ThermoParcel<ParcelType>::calc
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -138,6 +138,7 @@ Foam::layerParameters::layerParameters
readLabel(dict.lookup("nSmoothSurfaceNormals"))
),
nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
nSmoothDisplacement_(dict.lookupOrDefault("nSmoothDisplacement", 0)),
nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
maxFaceThicknessRatio_
(
@ -278,7 +279,7 @@ Foam::layerParameters::layerParameters
const keyType& key = iter().keyword();
const labelHashSet patchIDs
(
boundaryMesh.patchSet(List<wordRe>(1, key))
boundaryMesh.patchSet(List<wordRe>(1, wordRe(key)))
);
if (patchIDs.size() == 0)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -119,6 +119,8 @@ private:
label nSmoothNormals_;
label nSmoothDisplacement_;
label nSmoothThickness_;
scalar maxFaceThicknessRatio_;
@ -275,6 +277,12 @@ public:
return layerTerminationCos_;
}
//- Smooth internal displacement
label nSmoothDisplacement() const
{
return nSmoothDisplacement_;
}
//- Smooth layer thickness over surface patches
label nSmoothThickness() const
{

View File

@ -50,14 +50,11 @@ SourceFiles
#define AMIInterpolation_H
#include "className.H"
#include "DynamicList.H"
#include "searchableSurface.H"
#include "treeBoundBoxList.H"
#include "boolList.H"
#include "primitivePatch.H"
#include "faceAreaIntersect.H"
#include "indexedOctree.H"
#include "treeDataPrimitivePatch.H"
#include "treeBoundBoxList.H"
#include "globalIndex.H"
#include "ops.H"
@ -82,12 +79,38 @@ class AMIInterpolation
:
public AMIInterpolationName
{
//- local typedef to octree tree-type
typedef treeDataPrimitivePatch<TargetPatch> treeType;
public:
// Public data types
//- Enumeration specifying interpolation method
enum interpolationMethod
{
imDirect,
imMapNearest,
imFaceAreaWeight
};
//- Convert interpolationMethod to word representation
static word interpolationMethodToWord
(
const interpolationMethod& method
);
//- Convert word to interpolationMethod
static interpolationMethod wordTointerpolationMethod
(
const word& method
);
private:
// Private data
//- Interpolation method
interpolationMethod method_;
//- Flag to indicate that the two patches are co-directional and
// that the orientation of the target patch should be reversed
const bool reverseTarget_;
@ -96,6 +119,7 @@ class AMIInterpolation
// cases
label singlePatchProc_;
// Source patch
//- Source face areas
@ -110,10 +134,6 @@ class AMIInterpolation
//- Sum of weights of target faces per source face
scalarField srcWeightsSum_;
//- Labels of faces that are not overlapped by any target faces
// (should be empty for correct functioning)
labelList srcNonOverlap_;
// Target patch
@ -130,9 +150,6 @@ class AMIInterpolation
scalarField tgtWeightsSum_;
//- Octree used to find face seeds
autoPtr<indexedOctree<treeType> > treePtr_;
//- Face triangulation mode
const faceAreaIntersect::triangulationMode triMode_;
@ -152,30 +169,6 @@ class AMIInterpolation
void operator=(const AMIInterpolation&);
// Helper functions
//- Write triangle intersection to OBJ file
void writeIntersectionOBJ
(
const scalar area,
const face& f1,
const face& f2,
const pointField& f1Points,
const pointField& f2Points
) const;
//- Check that patches are valid
void checkPatches
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const;
//- Reset the octree for the traget patch face search
void resetTree(const TargetPatch& tgtPatch);
// Parallel functionality
//- Calculate if patches are on multiple processors
@ -229,83 +222,8 @@ class AMIInterpolation
) const;
// Marching front
//- Find face on target patch that overlaps source face
label findTargetFace
(
const label srcFaceI,
const SourcePatch& srcPatch
) const;
//- Add faces neighbouring faceI to the ID list
void appendNbrFaces
(
const label faceI,
const TargetPatch& patch,
const DynamicList<label>& visitedFaces,
DynamicList<label>& faceIDs
) const;
bool processSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const label srcFaceI,
const label tgtStartFaceI,
DynamicList<label>& nbrFaces,
DynamicList<label>& visitedFaces,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
void restartUncoveredSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
//- Set the source and target seed faces
void setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const SourcePatch& srcPatch0,
const TargetPatch& tgtPatch0,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI,
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const;
//- Calculate addressing
void calcAddressing
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
label srcFaceI = -1,
label tgtFaceI = -1
);
//- Normalise the (area) weights - suppresses numerical error in
// weights calculation
// NOTE: if area weights are incorrect by 'a significant amount'
@ -352,6 +270,7 @@ public:
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const faceAreaIntersect::triangulationMode& triMode,
const interpolationMethod& method = imFaceAreaWeight,
const bool reverseTarget = false
);
@ -362,6 +281,7 @@ public:
const TargetPatch& tgtPatch,
const autoPtr<searchableSurface>& surf,
const faceAreaIntersect::triangulationMode& triMode,
const interpolationMethod& method = imFaceAreaWeight,
const bool reverseTarget = false
);
@ -378,6 +298,12 @@ public:
//- Destructor
~AMIInterpolation();
// Typedef to SourcePatch type this AMIInterplation is instantiated on
typedef SourcePatch sourcePatchType;
// Typedef to TargetPatch type this AMIInterplation is instantiated on
typedef TargetPatch targetPatchType;
// Member Functions
@ -403,10 +329,6 @@ public:
// patch weights (i.e. the sum before normalisation)
inline const scalarField& srcWeightsSum() const;
//- Labels of faces that are not overlapped by any target faces
// (should be empty for correct functioning)
inline const labelList& srcNonOverlap() const;
//- Source map pointer - valid only if singlePatchProc = -1
// This gets source data into a form to be consumed by
// tgtAddress, tgtWeights

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -55,14 +55,6 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeights() const
}
template<class SourcePatch, class TargetPatch>
inline const Foam::labelList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcNonOverlap() const
{
return srcNonOverlap_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarField&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum() const

View File

@ -0,0 +1,339 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "AMIMethod.H"
#include "meshTools.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::AMIMethod<SourcePatch, TargetPatch>::checkPatches() const
{
if (debug && (!srcPatch_.size() || !tgtPatch_.size()))
{
Pout<< "AMI: Patches not on processor: Source faces = "
<< srcPatch_.size() << ", target faces = " << tgtPatch_.size()
<< endl;
}
const scalar maxBoundsError = 0.05;
// check bounds of source and target
boundBox bbSrc(srcPatch_.points(), srcPatch_.meshPoints(), true);
boundBox bbTgt(tgtPatch_.points(), tgtPatch_.meshPoints(), true);
boundBox bbTgtInf(bbTgt);
bbTgtInf.inflate(maxBoundsError);
if (!bbTgtInf.contains(bbSrc))
{
WarningIn("AMIMethod<SourcePatch, TargetPatch>::checkPatches()")
<< "Source and target patch bounding boxes are not similar" << nl
<< " source box span : " << bbSrc.span() << nl
<< " target box span : " << bbTgt.span() << nl
<< " source box : " << bbSrc << nl
<< " target box : " << bbTgt << nl
<< " inflated target box : " << bbTgtInf << endl;
}
}
template<class SourcePatch, class TargetPatch>
bool Foam::AMIMethod<SourcePatch, TargetPatch>::initialise
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label& srcFaceI,
label& tgtFaceI
)
{
// set initial sizes for weights and addressing - must be done even if
// returns false below
srcAddress.setSize(srcPatch_.size());
srcWeights.setSize(srcPatch_.size());
tgtAddress.setSize(tgtPatch_.size());
tgtWeights.setSize(tgtPatch_.size());
// check that patch sizes are valid
if (!srcPatch_.size())
{
return false;
}
else if (!tgtPatch_.size())
{
WarningIn
(
"void Foam::AMIMethod<SourcePatch, TargetPatch>::initialise"
"("
"label&, "
"label&"
")"
)
<< srcPatch_.size() << " source faces but no target faces" << endl;
return false;
}
// reset the octree
resetTree();
// find initial face match using brute force/octree search
if ((srcFaceI == -1) || (tgtFaceI == -1))
{
srcFaceI = 0;
tgtFaceI = 0;
bool foundFace = false;
forAll(srcPatch_, faceI)
{
tgtFaceI = findTargetFace(faceI);
if (tgtFaceI >= 0)
{
srcFaceI = faceI;
foundFace = true;
break;
}
}
if (!foundFace)
{
FatalErrorIn
(
"void Foam::AMIMethod<SourcePatch, TargetPatch>::initialise"
"("
"label&, "
"label&"
")"
) << "Unable to find initial target face" << abort(FatalError);
}
}
if (debug)
{
Pout<< "AMI: initial target face = " << tgtFaceI << endl;
}
return true;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIMethod<SourcePatch, TargetPatch>::writeIntersectionOBJ
(
const scalar area,
const face& f1,
const face& f2,
const pointField& f1Points,
const pointField& f2Points
) const
{
static label count = 1;
const pointField f1pts = f1.points(f1Points);
const pointField f2pts = f2.points(f2Points);
Pout<< "Face intersection area (" << count << "):" << nl
<< " f1 face = " << f1 << nl
<< " f1 pts = " << f1pts << nl
<< " f2 face = " << f2 << nl
<< " f2 pts = " << f2pts << nl
<< " area = " << area
<< endl;
OFstream os("areas" + name(count) + ".obj");
forAll(f1pts, i)
{
meshTools::writeOBJ(os, f1pts[i]);
}
os<< "l";
forAll(f1pts, i)
{
os<< " " << i + 1;
}
os<< " 1" << endl;
forAll(f2pts, i)
{
meshTools::writeOBJ(os, f2pts[i]);
}
os<< "l";
forAll(f2pts, i)
{
os<< " " << f1pts.size() + i + 1;
}
os<< " " << f1pts.size() + 1 << endl;
count++;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree()
{
// Clear the old octree
treePtr_.clear();
treeBoundBox bb(tgtPatch_.points());
bb.inflate(0.01);
if (!treePtr_.valid())
{
treePtr_.reset
(
new indexedOctree<treeType>
(
treeType(false, tgtPatch_),
bb, // overall search domain
8, // maxLevel
10, // leaf size
3.0 // duplicity
)
);
}
}
template<class SourcePatch, class TargetPatch>
Foam::label Foam::AMIMethod<SourcePatch, TargetPatch>::findTargetFace
(
const label srcFaceI
) const
{
label targetFaceI = -1;
const pointField& srcPts = srcPatch_.points();
const face& srcFace = srcPatch_[srcFaceI];
const point srcPt = srcFace.centre(srcPts);
const scalar srcFaceArea = srcMagSf_[srcFaceI];
pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
if (debug)
{
Pout<< "Source point = " << srcPt << ", Sample point = "
<< sample.hitPoint() << ", Sample index = " << sample.index()
<< endl;
}
if (sample.hit())
{
targetFaceI = sample.index();
}
return targetFaceI;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIMethod<SourcePatch, TargetPatch>::appendNbrFaces
(
const label faceI,
const TargetPatch& patch,
const DynamicList<label>& visitedFaces,
DynamicList<label>& faceIDs
) const
{
const labelList& nbrFaces = patch.faceFaces()[faceI];
// filter out faces already visited from src face neighbours
forAll(nbrFaces, i)
{
label nbrFaceI = nbrFaces[i];
bool valid = true;
forAll(visitedFaces, j)
{
if (nbrFaceI == visitedFaces[j])
{
valid = false;
break;
}
}
if (valid)
{
forAll(faceIDs, j)
{
if (nbrFaceI == faceIDs[j])
{
valid = false;
break;
}
}
}
if (valid)
{
faceIDs.append(nbrFaceI);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::AMIMethod<SourcePatch, TargetPatch>::AMIMethod
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
srcPatch_(srcPatch),
tgtPatch_(tgtPatch),
reverseTarget_(reverseTarget),
srcMagSf_(srcMagSf),
tgtMagSf_(tgtMagSf),
srcNonOverlap_(),
triMode_(triMode)
{
checkPatches();
label srcSize = returnReduce(srcPatch_.size(), sumOp<label>());
label tgtSize = returnReduce(tgtPatch_.size(), sumOp<label>());
IInfo<< "AMI: Creating addressing and weights between "
<< srcSize << " source faces and " << tgtSize << " target faces"
<< endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::AMIMethod<SourcePatch, TargetPatch>::~AMIMethod()
{}
// ************************************************************************* //

View File

@ -0,0 +1,268 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::AMIMethod
Description
Base class for Arbitrary Mesh Interface (AMI) methods
SourceFiles
AMIMethod.C
\*---------------------------------------------------------------------------*/
#ifndef AMIMethod_H
#define AMIMethod_H
#include "className.H"
#include "DynamicList.H"
#include "faceAreaIntersect.H"
#include "indexedOctree.H"
#include "treeDataPrimitivePatch.H"
#include "treeBoundBoxList.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class AMIMethod Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class AMIMethod
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
AMIMethod(const AMIMethod&);
//- Disallow default bitwise assignment
void operator=(const AMIMethod&);
protected:
//- local typedef to octree tree-type
typedef treeDataPrimitivePatch<TargetPatch> treeType;
// Protected data
//- Reference to source patch
const SourcePatch& srcPatch_;
//- Reference to target patch
const SourcePatch& tgtPatch_;
//- Flag to indicate that the two patches are co-directional and
// that the orientation of the target patch should be reversed
const bool reverseTarget_;
//- Source face areas
const scalarField& srcMagSf_;
//- Target face areas
const scalarField& tgtMagSf_;
//- Labels of faces that are not overlapped by any target faces
// (should be empty for correct functioning)
labelList srcNonOverlap_;
//- Octree used to find face seeds
autoPtr<indexedOctree<treeType> > treePtr_;
//- Face triangulation mode
const faceAreaIntersect::triangulationMode triMode_;
// Protected Member Functions
// Helper functions
//- Check AMI patch coupling
void checkPatches() const;
//- Initialise and return true if all ok
bool initialise
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label& srcFaceI,
label& tgtFaceI
);
//- Write triangle intersection to OBJ file
void writeIntersectionOBJ
(
const scalar area,
const face& f1,
const face& f2,
const pointField& f1Points,
const pointField& f2Points
) const;
// Common AMI method functions
//- Reset the octree for the target patch face search
void resetTree();
//- Find face on target patch that overlaps source face
label findTargetFace(const label srcFaceI) const;
//- Add faces neighbouring faceI to the ID list
void appendNbrFaces
(
const label faceI,
const TargetPatch& patch,
const DynamicList<label>& visitedFaces,
DynamicList<label>& faceIDs
) const;
public:
//- Runtime type information
TypeName("AMIMethod");
//- Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
AMIMethod,
components,
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
),
(srcPatch, tgtPatch, srcMagSf, tgtMagSf, triMode, reverseTarget)
);
//- Construct from components
AMIMethod
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
);
//- Selector
static autoPtr<AMIMethod> New
(
const word& methodName,
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
);
//- Destructor
virtual ~AMIMethod();
// Member Functions
// Access
//- Labels of faces that are not overlapped by any target faces
// Note: this should be empty for correct functioning
inline const labelList& srcNonOverlap() const;
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
) = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeAMIMethod(AMIType) \
\
typedef AMIMethod<AMIType::sourcePatchType,AMIType::targetPatchType> \
AMIMethod##AMIType; \
\
defineNamedTemplateTypeNameAndDebug(AMIMethod##AMIType, 0); \
defineTemplateRunTimeSelectionTable(AMIMethod##AMIType, components);
#define makeAMIMethodType(AMIType, Method) \
\
typedef Method<AMIType::sourcePatchType,AMIType::targetPatchType> \
Method##AMIType; \
\
defineNamedTemplateTypeNameAndDebug(Method##AMIType, 0); \
\
AMIMethod<AMIType::sourcePatchType,AMIType::targetPatchType>:: \
addcomponentsConstructorToTable<Method##AMIType> \
add##Method##AMIType##ConstructorToTable_;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "AMIMethodI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "AMIMethod.C"
# include "AMIMethodNew.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,34 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
inline const Foam::labelList&
Foam::AMIMethod<SourcePatch, TargetPatch>::srcNonOverlap() const
{
return srcNonOverlap_;
}
// ************************************************************************* //

View File

@ -0,0 +1,84 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::autoPtr<Foam::AMIMethod<SourcePatch, TargetPatch> >
Foam::AMIMethod<SourcePatch, TargetPatch>::New
(
const word& methodName,
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
{
if (debug)
{
Info<< "Selecting AMIMethod " << methodName << endl;
}
typename componentsConstructorTable::iterator cstrIter =
componentsConstructorTablePtr_->find(methodName);
if (cstrIter == componentsConstructorTablePtr_->end())
{
FatalErrorIn
(
"AMIMethod<SourcePatch, TargetPatch>::New"
"("
"const word&, "
"const SourcePatch&, "
"const TargetPatch&, "
"const scalarField&, "
"const scalarField&, "
"const faceAreaIntersect::triangulationMode&, "
"const bool"
")"
) << "Unknown AMIMethod type "
<< methodName << nl << nl
<< "Valid AMIMethod types are:" << nl
<< componentsConstructorTablePtr_->sortedToc() << exit(FatalError);
}
return autoPtr<AMIMethod<SourcePatch, TargetPatch> >
(
cstrIter()
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
);
}
// ************************************************************************* //

View File

@ -0,0 +1,224 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "directAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
(
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
DynamicList<label>& nonOverlapFaces,
label& srcFaceI,
label& tgtFaceI
) const
{
const labelList& srcNbr = this->srcPatch_.faceFaces()[srcFaceI];
const labelList& tgtNbr = this->tgtPatch_.faceFaces()[tgtFaceI];
const pointField& srcPoints = this->srcPatch_.points();
const pointField& tgtPoints = this->tgtPatch_.points();
const vectorField& srcCf = this->srcPatch_.faceCentres();
forAll(srcNbr, i)
{
label srcI = srcNbr[i];
if (mapFlag[srcI] && srcTgtSeed[srcI] == -1)
{
const face& srcF = this->srcPatch_[srcI];
const vector srcN = srcF.normal(srcPoints);
bool found = false;
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
const face& tgtF = this->tgtPatch_[tgtI];
pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
if (ray.hit())
{
// new match - append to lists
found = true;
srcTgtSeed[srcI] = tgtI;
srcSeeds.append(srcI);
break;
}
}
if (!found)
{
// no match available for source face srcI
mapFlag[srcI] = false;
nonOverlapFaces.append(srcI);
}
}
}
if (srcSeeds.size())
{
srcFaceI = srcSeeds.remove();
tgtFaceI = srcTgtSeed[srcFaceI];
}
else
{
srcFaceI = -1;
tgtFaceI = -1;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::directAMI<SourcePatch, TargetPatch>::directAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
AMIMethod<SourcePatch, TargetPatch>
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::directAMI<SourcePatch, TargetPatch>::~directAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::directAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI,
label tgtFaceI
)
{
bool ok =
this->initialise
(
srcAddress,
srcWeights,
tgtAddress,
tgtWeights,
srcFaceI,
tgtFaceI
);
if (!ok)
{
return;
}
// temporary storage for addressing and weights
List<DynamicList<label> > srcAddr(this->srcPatch_.size());
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> srcSeeds(10);
// list to keep track of tgt faces used to seed src faces
labelList srcTgtSeed(srcAddr.size(), -1);
srcTgtSeed[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(srcAddr.size(), true);
DynamicList<label> nonOverlapFaces;
do
{
srcAddr[srcFaceI].append(tgtFaceI);
tgtAddr[tgtFaceI].append(srcFaceI);
mapFlag[srcFaceI] = false;
// Do advancing front starting from srcFaceI, tgtFaceI
appendToDirectSeeds
(
mapFlag,
srcTgtSeed,
srcSeeds,
nonOverlapFaces,
srcFaceI,
tgtFaceI
);
} while (srcFaceI >= 0);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
<< " non-overlap faces identified"
<< endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
// transfer data to persistent storage
forAll(srcAddr, i)
{
scalar magSf = this->srcMagSf_[i];
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i] = scalarList(1, magSf);
}
forAll(tgtAddr, i)
{
scalar magSf = this->tgtMagSf_[i];
tgtAddress[i].transfer(tgtAddr[i]);
tgtWeights[i] = scalarList(1, magSf);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,144 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::directAMI
Description
Direct mapped Arbitrary Mesh Interface (AMI) method
SourceFiles
directAMI.C
\*---------------------------------------------------------------------------*/
#ifndef directAMI_H
#define directAMI_H
#include "AMIMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class directAMI Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class directAMI
:
public AMIMethod<SourcePatch, TargetPatch>
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
directAMI(const directAMI&);
//- Disallow default bitwise assignment
void operator=(const directAMI&);
// Marching front
//- Append to list of src face seed indices
void appendToDirectSeeds
(
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
DynamicList<label>& nonOverlapFaces,
label& srcFaceI,
label& tgtFaceI
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI
) const;
public:
//- Runtime type information
TypeName("directAMI");
// Constructors
//- Construct from components
directAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget = false
);
//- Destructor
virtual ~directAMI();
// Member Functions
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "directAMI.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,553 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
(
const label srcFaceI,
const label tgtStartFaceI,
// list of tgt face neighbour faces
DynamicList<label>& nbrFaces,
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label>& visitedFaces,
// temporary storage for addressing and weights
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
)
{
nbrFaces.clear();
visitedFaces.clear();
// append initial target face and neighbours
nbrFaces.append(tgtStartFaceI);
this->appendNbrFaces
(
tgtStartFaceI,
this->tgtPatch_,
visitedFaces,
nbrFaces
);
bool faceProcessed = false;
do
{
// process new target face
label tgtFaceI = nbrFaces.remove();
visitedFaces.append(tgtFaceI);
scalar area = interArea(srcFaceI, tgtFaceI);
// store when intersection area > 0
if (area > 0)
{
srcAddr[srcFaceI].append(tgtFaceI);
srcWght[srcFaceI].append(area);
tgtAddr[tgtFaceI].append(srcFaceI);
tgtWght[tgtFaceI].append(area);
this->appendNbrFaces
(
tgtFaceI,
this->tgtPatch_,
visitedFaces,
nbrFaces
);
faceProcessed = true;
}
} while (nbrFaces.size() > 0);
return faceProcessed;
}
template<class SourcePatch, class TargetPatch>
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
) const
{
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI];
// set possible seeds for later use
bool valuesSet = false;
forAll(srcNbrFaces, i)
{
label faceS = srcNbrFaces[i];
if (mapFlag[faceS] && seedFaces[faceS] == -1)
{
forAll(visitedFaces, j)
{
label faceT = visitedFaces[j];
scalar area = interArea(faceS, faceT);
scalar areaTotal = this->srcMagSf_[srcFaceI];
// Check that faces have enough overlap for robust walking
if (area/areaTotal > faceAreaIntersect::tolerance())
{
// TODO - throwing area away - re-use in next iteration?
seedFaces[faceS] = faceT;
if (!valuesSet)
{
srcFaceI = faceS;
tgtFaceI = faceT;
valuesSet = true;
}
}
}
}
}
// set next src and tgt faces if not set above
if (valuesSet)
{
return;
}
else
{
// try to use existing seed
bool foundNextSeed = false;
for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
{
if (mapFlag[faceI])
{
if (!foundNextSeed)
{
startSeedI = faceI;
foundNextSeed = true;
}
if (seedFaces[faceI] != -1)
{
srcFaceI = faceI;
tgtFaceI = seedFaces[faceI];
return;
}
}
}
// perform new search to find match
if (debug)
{
Pout<< "Advancing front stalled: searching for new "
<< "target face" << endl;
}
foundNextSeed = false;
for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
{
if (mapFlag[faceI])
{
if (!foundNextSeed)
{
startSeedI = faceI + 1;
foundNextSeed = true;
}
srcFaceI = faceI;
tgtFaceI = this->findTargetFace(srcFaceI);
if (tgtFaceI >= 0)
{
return;
}
}
}
FatalErrorIn
(
"void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::"
"setNextFaces"
"("
"label&, "
"label&, "
"label&, "
"const boolList&, "
"labelList&, "
"const DynamicList<label>&"
") const"
) << "Unable to set source and target faces" << abort(FatalError);
}
}
template<class SourcePatch, class TargetPatch>
Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
(
const label srcFaceI,
const label tgtFaceI
) const
{
const pointField& srcPoints = this->srcPatch_.points();
const pointField& tgtPoints = this->tgtPatch_.points();
// references to candidate faces
const face& src = this->srcPatch_[srcFaceI];
const face& tgt = this->tgtPatch_[tgtFaceI];
// quick reject if either face has zero area
// Note: do not used stored face areas for target patch
if
(
(this->srcMagSf_[srcFaceI] < ROOTVSMALL)
|| (tgt.mag(tgtPoints) < ROOTVSMALL)
)
{
return 0.0;
}
// create intersection object
faceAreaIntersect inter(srcPoints, tgtPoints, this->reverseTarget_);
// crude resultant norm
vector n(-src.normal(srcPoints));
if (this->reverseTarget_)
{
n -= tgt.normal(tgtPoints);
}
else
{
n += tgt.normal(tgtPoints);
}
n *= 0.5;
scalar area = 0;
if (mag(n) > ROOTVSMALL)
{
area = inter.calc(src, tgt, n, this->triMode_);
}
else
{
WarningIn
(
"void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::"
"interArea"
"("
"const label, "
"const label"
") const"
) << "Invalid normal for source face " << srcFaceI
<< " points " << UIndirectList<point>(srcPoints, src)
<< " target face " << tgtFaceI
<< " points " << UIndirectList<point>(tgtPoints, tgt)
<< endl;
}
if ((debug > 1) && (area > 0))
{
this->writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
}
return area;
}
template<class SourcePatch, class TargetPatch>
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::
restartUncoveredSourceFace
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
)
{
// Collect all src faces with a low weight
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelHashSet lowWeightFaces(100);
forAll(srcWght, srcFaceI)
{
scalar s = sum(srcWght[srcFaceI]);
scalar t = s/this->srcMagSf_[srcFaceI];
if (t < 0.5)
{
lowWeightFaces.insert(srcFaceI);
}
}
if (debug)
{
Pout<< "faceAreaWeightAMI: restarting search on "
<< lowWeightFaces.size() << " faces since sum of weights < 0.5"
<< endl;
}
if (lowWeightFaces.size() > 0)
{
// Erase all the lowWeight source faces from the target
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DynamicList<label> okSrcFaces(10);
DynamicList<scalar> okSrcWeights(10);
forAll(tgtAddr, tgtFaceI)
{
okSrcFaces.clear();
okSrcWeights.clear();
DynamicList<label>& srcFaces = tgtAddr[tgtFaceI];
DynamicList<scalar>& srcWeights = tgtWght[tgtFaceI];
forAll(srcFaces, i)
{
if (!lowWeightFaces.found(srcFaces[i]))
{
okSrcFaces.append(srcFaces[i]);
okSrcWeights.append(srcWeights[i]);
}
}
if (okSrcFaces.size() < srcFaces.size())
{
srcFaces.transfer(okSrcFaces);
srcWeights.transfer(okSrcWeights);
}
}
// Restart search from best hit
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
forAllConstIter(labelHashSet, lowWeightFaces, iter)
{
label srcFaceI = iter.key();
label tgtFaceI = this->findTargetFace(srcFaceI);
if (tgtFaceI != -1)
{
//bool faceProcessed =
processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
// ? Check faceProcessed to see if restarting has worked.
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::faceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
AMIMethod<SourcePatch, TargetPatch>
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::~faceAreaWeightAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI,
label tgtFaceI
)
{
bool ok =
this->initialise
(
srcAddress,
srcWeights,
tgtAddress,
tgtWeights,
srcFaceI,
tgtFaceI
);
if (!ok)
{
return;
}
// temporary storage for addressing and weights
List<DynamicList<label> > srcAddr(this->srcPatch_.size());
List<DynamicList<scalar> > srcWght(srcAddr.size());
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size());
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcAddr.size();
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
<< " non-overlap faces identified"
<< endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
// Check for badly covered faces
if (debug)
{
restartUncoveredSourceFace
(
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
}
// transfer data to persistent storage
forAll(srcAddr, i)
{
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i].transfer(srcWght[i]);
}
forAll(tgtAddr, i)
{
tgtAddress[i].transfer(tgtAddr[i]);
tgtWeights[i].transfer(tgtWght[i]);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,166 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::faceAreaWeightAMI
Description
Face area weighted Arbitrary Mesh Interface (AMI) method
SourceFiles
faceAreaWeightAMI.C
\*---------------------------------------------------------------------------*/
#ifndef faceAreaWeightAMI_H
#define faceAreaWeightAMI_H
#include "AMIMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class faceAreaWeightAMI Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class faceAreaWeightAMI
:
public AMIMethod<SourcePatch, TargetPatch>
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
faceAreaWeightAMI(const faceAreaWeightAMI&);
//- Disallow default bitwise assignment
void operator=(const faceAreaWeightAMI&);
// Marching front
//- Determine overlap contributions for source face srcFaceI
bool processSourceFace
(
const label srcFaceI,
const label tgtStartFaceI,
DynamicList<label>& nbrFaces,
DynamicList<label>& visitedFaces,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
//- Attempt to re-evaluate source faces that have not been included
void restartUncoveredSourceFace
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
//- Set the source and target seed faces
void setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI
) const;
public:
//- Runtime type information
TypeName("faceAreaWeightAMI");
// Constructors
//- Construct from components
faceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget = false
);
//- Destructor
virtual ~faceAreaWeightAMI();
// Member Functions
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "faceAreaWeightAMI.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,343 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "mapNearestAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::mapNearestAMI<SourcePatch, TargetPatch>::findNearestFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const label& srcFaceI,
label& tgtFaceI
) const
{
const vectorField& srcCf = srcPatch.faceCentres();
const vectorField& tgtCf = tgtPatch.faceCentres();
const vector srcP = srcCf[srcFaceI];
DynamicList<label> tgtFaces(10);
tgtFaces.append(tgtFaceI);
DynamicList<label> visitedFaces(10);
scalar d = GREAT;
do
{
label tgtI = tgtFaces.remove();
visitedFaces.append(tgtI);
scalar dTest = magSqr(tgtCf[tgtI] - srcP);
if (dTest < d)
{
tgtFaceI = tgtI;
d = dTest;
this->appendNbrFaces
(
tgtFaceI,
tgtPatch,
visitedFaces,
tgtFaces
);
}
} while (tgtFaces.size() > 0);
}
template<class SourcePatch, class TargetPatch>
void Foam::mapNearestAMI<SourcePatch, TargetPatch>::setNextNearestFaces
(
boolList& mapFlag,
label& startSeedI,
label& srcFaceI,
label& tgtFaceI
) const
{
const labelList& srcNbr = this->srcPatch_.faceFaces()[srcFaceI];
srcFaceI = -1;
forAll(srcNbr, i)
{
label faceI = srcNbr[i];
if (mapFlag[faceI])
{
srcFaceI = faceI;
startSeedI = faceI + 1;
return;
}
}
forAll(mapFlag, srcFaceI)
{
if (!mapFlag[srcFaceI])
{
tgtFaceI = this->findTargetFace(srcFaceI);
if (tgtFaceI == -1)
{
const vectorField& srcCf = this->srcPatch_.faceCentres();
FatalErrorIn
(
"void Foam::mapNearestAMI<SourcePatch, TargetPatch>::"
"setNextNearestFaces"
"("
"boolList&, "
"label&, "
"label&, "
"label&"
") const"
)
<< "Unable to find target face for source face "
<< srcFaceI << " with face centre " << srcCf[srcFaceI]
<< abort(FatalError);
}
break;
}
}
}
template<class SourcePatch, class TargetPatch>
Foam::label Foam::mapNearestAMI<SourcePatch, TargetPatch>::findMappedSrcFace
(
const label tgtFaceI,
const List<DynamicList<label> >& tgtToSrc
) const
{
DynamicList<label> testFaces(10);
DynamicList<label> visitedFaces(10);
testFaces.append(tgtFaceI);
do
{
// search target tgtFaceI neighbours for match with source face
label tgtI = testFaces.remove();
if (findIndex(visitedFaces, tgtI) == -1)
{
visitedFaces.append(tgtI);
if (tgtToSrc[tgtI].size())
{
return tgtToSrc[tgtI][0];
}
else
{
const labelList& nbrFaces = this->tgtPatch_.faceFaces()[tgtI];
forAll(nbrFaces, i)
{
if (findIndex(visitedFaces, nbrFaces[i]) == -1)
{
testFaces.append(nbrFaces[i]);
}
}
}
}
} while (testFaces.size());
// did not find any match - should not be possible to get here!
return -1;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
AMIMethod<SourcePatch, TargetPatch>
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::mapNearestAMI<SourcePatch, TargetPatch>::~mapNearestAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::mapNearestAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI,
label tgtFaceI
)
{
bool ok =
this->initialise
(
srcAddress,
srcWeights,
tgtAddress,
tgtWeights,
srcFaceI,
tgtFaceI
);
if (!ok)
{
return;
}
// temporary storage for addressing and weights
List<DynamicList<label> > srcAddr(this->srcPatch_.size());
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// list to keep track of whether src face can be mapped
boolList mapFlag(srcAddr.size(), true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{
findNearestFace(this->srcPatch_, this->tgtPatch_, srcFaceI, tgtFaceI);
srcAddr[srcFaceI].append(tgtFaceI);
tgtAddr[tgtFaceI].append(srcFaceI);
mapFlag[srcFaceI] = false;
// Do advancing front starting from srcFaceI, tgtFaceI
setNextNearestFaces
(
mapFlag,
startSeedI,
srcFaceI,
tgtFaceI
);
} while (srcFaceI >= 0);
// for the case of multiple source faces per target face, select the
// nearest source face only and discard the others
const vectorField& srcCf = this->srcPatch_.faceCentres();
const vectorField& tgtCf = this->tgtPatch_.faceCentres();
forAll(tgtAddr, targetFaceI)
{
if (tgtAddr[targetFaceI].size() > 1)
{
const vector& tgtC = tgtCf[tgtFaceI];
DynamicList<label>& srcFaces = tgtAddr[targetFaceI];
label srcFaceI = srcFaces[0];
scalar d = magSqr(tgtC - srcCf[srcFaceI]);
for (label i = 1; i < srcFaces.size(); i++)
{
label srcI = srcFaces[i];
scalar dNew = magSqr(tgtC - srcCf[srcI]);
if (dNew < d)
{
d = dNew;
srcFaceI = srcI;
}
}
srcFaces.clear();
srcFaces.append(srcFaceI);
}
}
// If there are more target faces than source faces, some target faces
// might not yet be mapped
forAll(tgtAddr, tgtFaceI)
{
if (tgtAddr[tgtFaceI].empty())
{
label srcFaceI = findMappedSrcFace(tgtFaceI, tgtAddr);
// note - reversed search from src->tgt to tgt->src
findNearestFace
(
this->tgtPatch_,
this->srcPatch_,
tgtFaceI,
srcFaceI
);
tgtAddr[tgtFaceI].append(srcFaceI);
}
}
// transfer data to persistent storage
forAll(srcAddr, i)
{
scalar magSf = this->srcMagSf_[i];
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i] = scalarList(1, magSf);
}
forAll(tgtAddr, i)
{
scalar magSf = this->tgtMagSf_[i];
tgtAddress[i].transfer(tgtAddr[i]);
tgtWeights[i] = scalarList(1, magSf);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::mapNearestAMI
Description
Nearest-mapping Arbitrary Mesh Interface (AMI) method
SourceFiles
mapNearestAMI.C
\*---------------------------------------------------------------------------*/
#ifndef mapNearestAMI_H
#define mapNearestAMI_H
#include "AMIMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mapNearestAMI Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class mapNearestAMI
:
public AMIMethod<SourcePatch, TargetPatch>
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
mapNearestAMI(const mapNearestAMI&);
//- Disallow default bitwise assignment
void operator=(const mapNearestAMI&);
// Marching front
//- Find nearest target face for source face srcFaceI
void findNearestFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const label& srcFaceI,
label& tgtFaceI
) const;
//- Determine next source-target face pair
void setNextNearestFaces
(
boolList& mapFlag,
label& startSeedI,
label& srcFaceI,
label& tgtFaceI
) const;
//- Find mapped source face
label findMappedSrcFace
(
const label tgtFaceI,
const List<DynamicList<label> >& tgtToSrc
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI
) const;
public:
//- Runtime type information
TypeName("mapNearestAMI");
// Constructors
//- Construct from components
mapNearestAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget = false
);
//- Destructor
virtual ~mapNearestAMI();
// Member Functions
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "mapNearestAMI.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "AMIPatchToPatchInterpolation.H"
#include "AMIMethod.H"
#include "directAMI.H"
#include "mapNearestAMI.H"
#include "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makeAMIMethod(AMIPatchToPatchInterpolation);
makeAMIMethodType(AMIPatchToPatchInterpolation, directAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, mapNearestAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, faceAreaWeightAMI);
}
// ************************************************************************* //

View File

@ -283,6 +283,7 @@ void Foam::cyclicAMIPolyPatch::resetAMI() const
nbrPatch0,
surfPtr(),
faceAreaIntersect::tmMesh,
AMIPatchToPatchInterpolation::imFaceAreaWeight,
AMIReverse_
)
);

View File

@ -170,6 +170,7 @@ twoDPointCorrector/twoDPointCorrector.C
AMI=AMIInterpolation
$(AMI)/AMIInterpolation/AMIInterpolationName.C
$(AMI)/AMIInterpolation/AMIPatchToPatchInterpolation.C
$(AMI)/faceAreaIntersect/faceAreaIntersect.C
$(AMI)/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
$(AMI)/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -160,7 +160,11 @@ void Foam::treeDataPoint::findNearest
) const
{
// Best so far
scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
scalar nearestDistSqr = GREAT;
if (minIndex >= 0)
{
nearestDistSqr = magSqr(linePoint - nearestPoint);
}
forAll(indices, i)
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -404,11 +404,50 @@ void Foam::treeDataTriSurface::findNearest
point& nearestPoint
) const
{
notImplemented
(
"treeDataTriSurface::findNearest(const labelUList&"
", const linePointRef&, treeBoundBox&, label&, point&, point&) const"
);
// Best so far
scalar nearestDistSqr = VGREAT;
if (minIndex >= 0)
{
nearestDistSqr = magSqr(linePoint - nearestPoint);
}
const pointField& points = surface_.points();
forAll(indices, i)
{
label index = indices[i];
const triSurface::FaceType& f = surface_[index];
triPointRef tri(f.tri(points));
pointHit lineInfo(point::max);
pointHit pHit = tri.nearestPoint(ln, lineInfo);
scalar distSqr = sqr(pHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
linePoint = lineInfo.rawPoint();
nearestPoint = pHit.rawPoint();
{
point& minPt = tightest.min();
minPt = min(ln.start(), ln.end());
minPt.x() -= pHit.distance();
minPt.y() -= pHit.distance();
minPt.z() -= pHit.distance();
}
{
point& maxPt = tightest.max();
maxPt = max(ln.start(), ln.end());
maxPt.x() += pHit.distance();
maxPt.y() += pHit.distance();
maxPt.z() += pHit.distance();
}
}
}
}

View File

@ -849,6 +849,7 @@ void Foam::mappedPatchBase::calcAMI() const
samplePolyPatch(), // nbrPatch0,
surfPtr(),
faceAreaIntersect::tmMesh,
AMIPatchToPatchInterpolation::imFaceAreaWeight,
AMIReverse_
)
);

View File

@ -91,6 +91,7 @@ void Foam::regionCoupledBase::resetAMI() const
nbrPatch0,
surfPtr(),
faceAreaIntersect::tmMesh,
AMIPatchToPatchInterpolation::imFaceAreaWeight,
AMIReverse_
)
);

View File

@ -28,6 +28,7 @@ License
#include "faceZoneSet.H"
#include "searchableSurface.H"
#include "syncTools.H"
#include "Time.H"
#include "addToRunTimeSelectionTable.H"
@ -69,7 +70,15 @@ Foam::searchableSurfaceToFaceZone::searchableSurfaceToFaceZone
searchableSurface::New
(
word(dict.lookup("surface")),
IOobject
(
dict.lookupOrDefault("name", mesh.objectRegistry::db().name()),
mesh.time().constant(),
"triSurface",
mesh.objectRegistry::db(),
IOobject::MUST_READ,
IOobject::NO_WRITE
),
dict
)
)

View File

@ -221,6 +221,7 @@ Foam::regionModels::regionModel::interRegionAMI
p,
nbrP,
faceAreaIntersect::tmMesh,
AMIPatchToPatchInterpolation::imFaceAreaWeight,
flip
)
);
@ -261,6 +262,7 @@ Foam::regionModels::regionModel::interRegionAMI
p,
nbrP,
faceAreaIntersect::tmMesh,
AMIPatchToPatchInterpolation::imFaceAreaWeight,
flip
)
);
@ -318,7 +320,7 @@ Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID
(
"Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID"
"("
"const regionModel& , "
"const regionModel&, "
"const label"
") const"
)

View File

@ -60,10 +60,14 @@ $(meshToMesh)/calculateMeshToMeshAddressing.C
$(meshToMesh)/calculateMeshToMeshWeights.C
meshToMeshNew = meshToMeshInterpolation/meshToMeshNew
$(meshToMeshNew)/calcDirect.C
$(meshToMeshNew)/calcMapNearest.C
$(meshToMeshNew)/calcCellVolumeWeight.C
$(meshToMeshNew)/meshToMeshNew.C
$(meshToMeshNew)/meshToMeshNewParallelOps.C
meshToMeshNewMethods = meshToMeshInterpolation/meshToMeshNew/calcMethod
$(meshToMeshNewMethods)/meshToMeshMethod/meshToMeshMethod.C
$(meshToMeshNewMethods)/meshToMeshMethod/meshToMeshMethodNew.C
$(meshToMeshNewMethods)/cellVolumeWeight/cellVolumeWeightMethod.C
$(meshToMeshNewMethods)/direct/directMethod.C
$(meshToMeshNewMethods)/mapNearest/mapNearestMethod.C
LIB = $(FOAM_LIBBIN)/libsampling

View File

@ -1,159 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-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 "meshToMeshNew.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::meshToMeshNew::calcDirect
(
const polyMesh& src,
const polyMesh& tgt,
const label srcSeedI,
const label tgtSeedI
)
{
// store a list of src cells already mapped
boolList srcSeedFlag(src.nCells(), true);
labelList srcTgtSeed(src.nCells(), -1);
List<DynamicList<label> > srcToTgt(src.nCells());
List<DynamicList<label> > tgtToSrc(tgt.nCells());
DynamicList<label> srcSeeds;
const scalarField& srcVc = src.cellVolumes();
const scalarField& tgtVc = tgt.cellVolumes();
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
do
{
// store src/tgt cell pair
srcToTgt[srcCellI].append(tgtCellI);
tgtToSrc[tgtCellI].append(srcCellI);
// mark source cell srcSeedI as matched
srcSeedFlag[srcCellI] = false;
// accumulate intersection volume
V_ += srcVc[srcCellI];
// find new source seed cell
appendToDirectSeeds
(
src,
tgt,
srcSeedFlag,
srcTgtSeed,
srcSeeds,
srcCellI,
tgtCellI
);
}
while (srcCellI >= 0);
// transfer addressing into persistent storage
forAll(srcToTgtCellAddr_, i)
{
scalar v = srcVc[i];
srcToTgtCellAddr_[i].transfer(srcToTgt[i]);
srcToTgtCellWght_[i] = scalarList(srcToTgtCellAddr_[i].size(), v);
}
forAll(tgtToSrcCellAddr_, i)
{
scalar v = tgtVc[i];
tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]);
tgtToSrcCellWght_[i] = scalarList(tgtToSrcCellAddr_[i].size(), v);
}
}
void Foam::meshToMeshNew::appendToDirectSeeds
(
const polyMesh& src,
const polyMesh& tgt,
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
label& srcSeedI,
label& tgtSeedI
) const
{
const labelList& srcNbr = src.cellCells()[srcSeedI];
const labelList& tgtNbr = tgt.cellCells()[tgtSeedI];
const vectorField& srcCentre = src.cellCentres();
forAll(srcNbr, i)
{
label srcI = srcNbr[i];
if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1))
{
// source cell srcI not yet mapped
// identfy if target cell exists for source cell srcI
bool found = false;
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
if (tgt.pointInCell(srcCentre[srcI], tgtI))
{
// new match - append to lists
found = true;
srcTgtSeed[srcI] = tgtI;
srcSeeds.append(srcI);
break;
}
}
if (!found)
{
// no match available for source cell srcI
mapFlag[srcI] = false;
}
}
}
if (srcSeeds.size())
{
srcSeedI = srcSeeds.remove();
tgtSeedI = srcTgtSeed[srcSeedI];
}
else
{
srcSeedI = -1;
tgtSeedI = -1;
}
}
// ************************************************************************* //

View File

@ -1,265 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-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 "meshToMeshNew.H"
#include "ListOps.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::meshToMeshNew::calcMapNearest
(
const polyMesh& src,
const polyMesh& tgt,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
)
{
List<DynamicList<label> > srcToTgt(src.nCells());
List<DynamicList<label> > tgtToSrc(tgt.nCells());
const scalarField& srcVc = src.cellVolumes();
const scalarField& tgtVc = tgt.cellVolumes();
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
do
{
// find nearest tgt cell
findNearestCell(src, tgt, srcCellI, tgtCellI);
// store src/tgt cell pair
srcToTgt[srcCellI].append(tgtCellI);
tgtToSrc[tgtCellI].append(srcCellI);
// mark source cell srcCellI and tgtCellI as matched
mapFlag[srcCellI] = false;
// accumulate intersection volume
V_ += srcVc[srcCellI];
// find new source cell
setNextNearestCells
(
startSeedI,
srcCellI,
tgtCellI,
mapFlag,
src,
tgt,
srcCellIDs
);
}
while (srcCellI >= 0);
// for the case of multiple source cells per target cell, select the
// nearest source cell only and discard the others
const vectorField& srcCc = src.cellCentres();
const vectorField& tgtCc = tgt.cellCentres();
forAll(tgtToSrc, targetCellI)
{
if (tgtToSrc[targetCellI].size() > 1)
{
const vector& tgtC = tgtCc[tgtCellI];
DynamicList<label>& srcCells = tgtToSrc[targetCellI];
label srcCellI = srcCells[0];
scalar d = magSqr(tgtC - srcCc[srcCellI]);
for (label i = 1; i < srcCells.size(); i++)
{
label srcI = srcCells[i];
scalar dNew = magSqr(tgtC - srcCc[srcI]);
if (dNew < d)
{
d = dNew;
srcCellI = srcI;
}
}
srcCells.clear();
srcCells.append(srcCellI);
}
}
// If there are more target cells than source cells, some target cells
// might not yet be mapped
forAll(tgtToSrc, tgtCellI)
{
if (tgtToSrc[tgtCellI].empty())
{
label srcCellI = findMappedSrcCell(tgt, tgtCellI, tgtToSrc);
findNearestCell(tgt, src, tgtCellI, srcCellI);
tgtToSrc[tgtCellI].append(srcCellI);
}
}
// transfer addressing into persistent storage
forAll(srcToTgtCellAddr_, i)
{
scalar v = srcVc[i];
srcToTgtCellWght_[i] = scalarList(srcToTgt[i].size(), v);
srcToTgtCellAddr_[i].transfer(srcToTgt[i]);
}
forAll(tgtToSrcCellAddr_, i)
{
scalar v = tgtVc[i];
tgtToSrcCellWght_[i] = scalarList(tgtToSrc[i].size(), v);
tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]);
}
}
void Foam::meshToMeshNew::findNearestCell
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
label& tgtCellI
)
{
const vectorField& srcC = src.cellCentres();
const vectorField& tgtC = tgt.cellCentres();
const vector& srcP = srcC[srcCellI];
DynamicList<label> tgtCells(10);
tgtCells.append(tgtCellI);
DynamicList<label> visitedCells(10);
scalar d = GREAT;
do
{
label tgtI = tgtCells.remove();
visitedCells.append(tgtI);
scalar dTest = magSqr(tgtC[tgtI] - srcP);
if (dTest < d)
{
tgtCellI = tgtI;
d = dTest;
appendNbrCells(tgtCellI, tgt, visitedCells, tgtCells);
}
} while (tgtCells.size() > 0);
}
void Foam::meshToMeshNew::setNextNearestCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
boolList& mapFlag,
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs
)
{
const labelList& srcNbr = src.cellCells()[srcCellI];
srcCellI = -1;
forAll(srcNbr, i)
{
label cellI = srcNbr[i];
if (mapFlag[cellI])
{
srcCellI = cellI;
startSeedI = cellI + 1;
return;
}
}
(void)findInitialSeeds
(
src,
tgt,
srcCellIDs,
mapFlag,
startSeedI,
srcCellI,
tgtCellI
);
}
Foam::label Foam::meshToMeshNew::findMappedSrcCell
(
const polyMesh& tgt,
const label tgtCellI,
const List<DynamicList<label> >& tgtToSrc
) const
{
DynamicList<label> testCells(10);
DynamicList<label> visitedCells(10);
testCells.append(tgtCellI);
do
{
// search target tgtCellI neighbours for match with source cell
label tgtI = testCells.remove();
if (findIndex(visitedCells, tgtI) == -1)
{
visitedCells.append(tgtI);
if (tgtToSrc[tgtI].size())
{
return tgtToSrc[tgtI][0];
}
else
{
const labelList& nbrCells = tgt.cellCells()[tgtI];
forAll(nbrCells, i)
{
if (findIndex(visitedCells, nbrCells[i]) == -1)
{
testCells.append(nbrCells[i]);
}
}
}
}
} while (testCells.size());
// did not find any match - should not be possible to get here!
return -1;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,15 +23,79 @@ License
\*---------------------------------------------------------------------------*/
#include "meshToMeshNew.H"
#include "tetOverlapVolume.H"
#include "cellVolumeWeightMethod.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
void Foam::meshToMeshNew::calcCellVolumeWeight
namespace Foam
{
defineTypeNameAndDebug(cellVolumeWeightMethod, 0);
addToRunTimeSelectionTable
(
meshToMeshMethod,
cellVolumeWeightMethod,
components
);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
bool Foam::cellVolumeWeightMethod::findInitialSeeds
(
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const
{
const cellList& srcCells = src_.cells();
const faceList& srcFaces = src_.faces();
const pointField& srcPts = src_.points();
for (label i = startSeedI; i < srcCellIDs.size(); i++)
{
label srcI = srcCellIDs[i];
if (mapFlag[srcI])
{
const pointField
pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
forAll(pts, ptI)
{
const point& pt = pts[ptI];
label tgtI = tgt_.cellTree().findInside(pt);
if (tgtI != -1 && intersect(srcI, tgtI))
{
srcSeedI = srcI;
tgtSeedI = tgtI;
return true;
}
}
}
}
if (debug)
{
Pout<< "could not find starting seed" << endl;
}
return false;
}
void Foam::cellVolumeWeightMethod::calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
@ -42,11 +106,11 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
List<DynamicList<label> > srcToTgtAddr(src.nCells());
List<DynamicList<scalar> > srcToTgtWght(src.nCells());
List<DynamicList<label> > srcToTgtAddr(src_.nCells());
List<DynamicList<scalar> > srcToTgtWght(src_.nCells());
List<DynamicList<label> > tgtToSrcAddr(tgt.nCells());
List<DynamicList<scalar> > tgtToSrcWght(tgt.nCells());
List<DynamicList<label> > tgtToSrcAddr(tgt_.nCells());
List<DynamicList<scalar> > tgtToSrcWght(tgt_.nCells());
// list of tgt cell neighbour cells
DynamicList<label> nbrTgtCells(10);
@ -55,10 +119,10 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
DynamicList<label> visitedTgtCells(10);
// list to keep track of tgt cells used to seed src cells
labelList seedCells(src.nCells(), -1);
labelList seedCells(src_.nCells(), -1);
seedCells[srcCellI] = tgtCellI;
const scalarField& srcVol = src.cellVolumes();
const scalarField& srcVol = src_.cellVolumes();
do
{
@ -67,14 +131,14 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
// append initial target cell and neighbours
nbrTgtCells.append(tgtCellI);
appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
appendNbrCells(tgtCellI, tgt_, visitedTgtCells, nbrTgtCells);
do
{
tgtCellI = nbrTgtCells.remove();
visitedTgtCells.append(tgtCellI);
scalar vol = interVol(src, tgt, srcCellI, tgtCellI);
scalar vol = interVol(srcCellI, tgtCellI);
// accumulate addressing and weights for valid intersection
if (vol/srcVol[srcCellI] > tolerance_)
@ -86,7 +150,7 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
tgtToSrcAddr[tgtCellI].append(srcCellI);
tgtToSrcWght[tgtCellI].append(vol);
appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
appendNbrCells(tgtCellI, tgt_, visitedTgtCells, nbrTgtCells);
// accumulate intersection volume
V_ += vol;
@ -102,8 +166,6 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
startSeedI,
srcCellI,
tgtCellI,
src,
tgt,
srcCellIDs,
mapFlag,
visitedTgtCells,
@ -113,34 +175,32 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
while (srcCellI != -1);
// transfer addressing into persistent storage
forAll(srcToTgtCellAddr_, i)
forAll(srcToTgtCellAddr, i)
{
srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]);
srcToTgtCellWght_[i].transfer(srcToTgtWght[i]);
srcToTgtCellAddr[i].transfer(srcToTgtAddr[i]);
srcToTgtCellWght[i].transfer(srcToTgtWght[i]);
}
forAll(tgtToSrcCellAddr_, i)
forAll(tgtToSrcCellAddr, i)
{
tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]);
tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]);
tgtToSrcCellAddr[i].transfer(tgtToSrcAddr[i]);
tgtToSrcCellWght[i].transfer(tgtToSrcWght[i]);
}
}
void Foam::meshToMeshNew::setNextCells
void Foam::cellVolumeWeightMethod::setNextCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs,
const boolList& mapFlag,
const DynamicList<label>& visitedCells,
labelList& seedCells
) const
{
const labelList& srcNbrCells = src.cellCells()[srcCellI];
const labelList& srcNbrCells = src_.cellCells()[srcCellI];
// set possible seeds for later use by querying all src cell neighbours
// with all visited target cells
@ -155,7 +215,7 @@ void Foam::meshToMeshNew::setNextCells
{
label cellT = visitedCells[j];
if (intersect(src, tgt, cellS, cellT))
if (intersect(cellS, cellT))
{
seedCells[cellS] = cellT;
@ -211,8 +271,6 @@ void Foam::meshToMeshNew::setNextCells
bool restart =
findInitialSeeds
(
src,
tgt,
srcCellIDs,
mapFlag,
startSeedI,
@ -233,68 +291,91 @@ void Foam::meshToMeshNew::setNextCells
}
bool Foam::meshToMeshNew::intersect
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cellVolumeWeightMethod::cellVolumeWeightMethod
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
const label tgtCellI
) const
{
scalar threshold = tolerance_*src.cellVolumes()[srcCellI];
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell
(
pointField
(
tgt.points(),
tgt.cellPoints()[tgtCellI]
)
);
return overlapEngine.cellCellOverlapMinDecomp
(
src,
srcCellI,
tgt,
tgtCellI,
bbTgtCell,
threshold
);
}
const polyMesh& tgt
)
:
meshToMeshMethod(src, tgt)
{}
Foam::scalar Foam::meshToMeshNew::interVol
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cellVolumeWeightMethod::~cellVolumeWeightMethod()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cellVolumeWeightMethod::calculate
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
const label tgtCellI
) const
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToSrcAddr,
scalarListList& tgtToSrcWght
)
{
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell
bool ok = initialise
(
pointField
(
tgt.points(),
tgt.cellPoints()[tgtCellI]
)
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght
);
scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
if (!ok)
{
return;
}
// (potentially) participating source mesh cells
const labelList srcCellIDs(maskCells());
// list to keep track of whether src cell can be mapped
boolList mapFlag(src_.nCells(), false);
UIndirectList<bool>(mapFlag, srcCellIDs) = true;
// find initial point in tgt mesh
label srcSeedI = -1;
label tgtSeedI = -1;
label startSeedI = 0;
bool startWalk =
findInitialSeeds
(
src,
srcCellI,
tgt,
tgtCellI,
bbTgtCell
srcCellIDs,
mapFlag,
startSeedI,
srcSeedI,
tgtSeedI
);
return vol;
if (startWalk)
{
calculateAddressing
(
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght,
srcSeedI,
tgtSeedI,
srcCellIDs,
mapFlag,
startSeedI
);
}
else
{
// if meshes are collocated, after inflating the source mesh bounding
// box tgt mesh cells may be transferred, but may still not overlap
// with the source mesh
return;
}
}

View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::cellVolumeWeightMethod
Description
Cell-volume-weighted mesh-to-mesh interpolation class
Volume conservative.
SourceFiles
cellVolumeWeightMethod.C
\*---------------------------------------------------------------------------*/
#ifndef cellVolumeWeightMethod_H
#define cellVolumeWeightMethod_H
#include "meshToMeshMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cellVolumeWeightMethod Declaration
\*---------------------------------------------------------------------------*/
class cellVolumeWeightMethod
:
public meshToMeshMethod
{
protected:
// Protected Member Functions
//- Find indices of overlapping cells in src and tgt meshes - returns
// true if found a matching pair
bool findInitialSeeds
(
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const;
//- Calculate the mesh-to-mesh addressing and weights
void calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
);
//- Set the next cells in the advancing front algorithm
void setNextCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
const labelList& srcCellIDs,
const boolList& mapFlag,
const DynamicList<label>& visitedCells,
labelList& seedCells
) const;
//- Disallow default bitwise copy construct
cellVolumeWeightMethod(const cellVolumeWeightMethod&);
//- Disallow default bitwise assignment
void operator=(const cellVolumeWeightMethod&);
public:
//- Run-time type information
TypeName("cellVolumeWeight");
//- Construct from source and target meshes
cellVolumeWeightMethod(const polyMesh& src, const polyMesh& tgt);
//- Destructor
virtual ~cellVolumeWeightMethod();
// Member Functions
// Evaluate
//- Calculate addressing and weights
virtual void calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToTgtAddr,
scalarListList& tgtToTgtWght
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,303 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "directMethod.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(directMethod, 0);
addToRunTimeSelectionTable(meshToMeshMethod, directMethod, components);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
bool Foam::directMethod::findInitialSeeds
(
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const
{
const cellList& srcCells = src_.cells();
const faceList& srcFaces = src_.faces();
const pointField& srcPts = src_.points();
for (label i = startSeedI; i < srcCellIDs.size(); i++)
{
label srcI = srcCellIDs[i];
if (mapFlag[srcI])
{
const pointField
pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
forAll(pts, ptI)
{
const point& pt = pts[ptI];
label tgtI = tgt_.cellTree().findInside(pt);
if (tgtI != -1 && intersect(srcI, tgtI))
{
srcSeedI = srcI;
tgtSeedI = tgtI;
return true;
}
}
}
}
if (debug)
{
Pout<< "could not find starting seed" << endl;
}
return false;
}
void Foam::directMethod::calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs, // not used
boolList& mapFlag,
label& startSeedI
)
{
// store a list of src cells already mapped
labelList srcTgtSeed(src_.nCells(), -1);
List<DynamicList<label> > srcToTgt(src_.nCells());
List<DynamicList<label> > tgtToSrc(tgt_.nCells());
DynamicList<label> srcSeeds(10);
const scalarField& srcVc = src_.cellVolumes();
const scalarField& tgtVc = tgt_.cellVolumes();
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
do
{
// store src/tgt cell pair
srcToTgt[srcCellI].append(tgtCellI);
tgtToSrc[tgtCellI].append(srcCellI);
// mark source cell srcSeedI as matched
mapFlag[srcCellI] = false;
// accumulate intersection volume
V_ += srcVc[srcCellI];
// find new source seed cell
appendToDirectSeeds
(
mapFlag,
srcTgtSeed,
srcSeeds,
srcCellI,
tgtCellI
);
}
while (srcCellI >= 0);
// transfer addressing into persistent storage
forAll(srcToTgtCellAddr, i)
{
scalar v = srcVc[i];
srcToTgtCellAddr[i].transfer(srcToTgt[i]);
srcToTgtCellWght[i] = scalarList(1, v);
}
forAll(tgtToSrcCellAddr, i)
{
scalar v = tgtVc[i];
tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
tgtToSrcCellWght[i] = scalarList(1, v);
}
}
void Foam::directMethod::appendToDirectSeeds
(
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
label& srcSeedI,
label& tgtSeedI
) const
{
const labelList& srcNbr = src_.cellCells()[srcSeedI];
const labelList& tgtNbr = tgt_.cellCells()[tgtSeedI];
const vectorField& srcCentre = src_.cellCentres();
forAll(srcNbr, i)
{
label srcI = srcNbr[i];
if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1))
{
// source cell srcI not yet mapped
// identfy if target cell exists for source cell srcI
bool found = false;
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
if (tgt_.pointInCell(srcCentre[srcI], tgtI))
{
// new match - append to lists
found = true;
srcTgtSeed[srcI] = tgtI;
srcSeeds.append(srcI);
break;
}
}
if (!found)
{
// no match available for source cell srcI
mapFlag[srcI] = false;
}
}
}
if (srcSeeds.size())
{
srcSeedI = srcSeeds.remove();
tgtSeedI = srcTgtSeed[srcSeedI];
}
else
{
srcSeedI = -1;
tgtSeedI = -1;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::directMethod::directMethod
(
const polyMesh& src,
const polyMesh& tgt
)
:
meshToMeshMethod(src, tgt)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::directMethod::~directMethod()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::directMethod::calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToSrcAddr,
scalarListList& tgtToSrcWght
)
{
bool ok = initialise
(
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght
);
if (!ok)
{
return;
}
// (potentially) participating source mesh cells
const labelList srcCellIDs(maskCells());
// list to keep track of whether src cell can be mapped
boolList mapFlag(src_.nCells(), false);
UIndirectList<bool>(mapFlag, srcCellIDs) = true;
// find initial point in tgt mesh
label srcSeedI = -1;
label tgtSeedI = -1;
label startSeedI = 0;
bool startWalk =
findInitialSeeds
(
srcCellIDs,
mapFlag,
startSeedI,
srcSeedI,
tgtSeedI
);
if (startWalk)
{
calculateAddressing
(
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght,
srcSeedI,
tgtSeedI,
srcCellIDs,
mapFlag,
startSeedI
);
}
else
{
// if meshes are collocated, after inflating the source mesh bounding
// box tgt mesh cells may be transferred, but may still not overlap
// with the source mesh
return;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::directMethod
Description
Direct (one-to-one cell correspondence) mesh-to-mesh interpolation class
Volume conservative.
SourceFiles
directMethod.C
\*---------------------------------------------------------------------------*/
#ifndef directMethod_H
#define directMethod_H
#include "meshToMeshMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class directMethod Declaration
\*---------------------------------------------------------------------------*/
class directMethod
:
public meshToMeshMethod
{
protected:
// Protected Member Functions
//- Find indices of overlapping cells in src and tgt meshes - returns
// true if found a matching pair
bool findInitialSeeds
(
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const;
//- Calculate the mesh-to-mesh addressing and weights
void calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
);
//- Append to list of src mesh seed indices
void appendToDirectSeeds
(
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
label& srcSeedI,
label& tgtSeedI
) const;
//- Disallow default bitwise copy construct
directMethod(const directMethod&);
//- Disallow default bitwise assignment
void operator=(const directMethod&);
public:
//- Run-time type information
TypeName("direct");
//- Construct from source and target meshes
directMethod(const polyMesh& src, const polyMesh& tgt);
//- Destructor
virtual ~directMethod();
// Member Functions
// Evaluate
//- Calculate addressing and weights
virtual void calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToTgtAddr,
scalarListList& tgtToTgtWght
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,424 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "mapNearestMethod.H"
#include "pointIndexHit.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(mapNearestMethod, 0);
addToRunTimeSelectionTable(meshToMeshMethod, mapNearestMethod, components);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
bool Foam::mapNearestMethod::findInitialSeeds
(
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const
{
const cellList& srcCells = src_.cells();
const faceList& srcFaces = src_.faces();
const pointField& srcPts = src_.points();
for (label i = startSeedI; i < srcCellIDs.size(); i++)
{
label srcI = srcCellIDs[i];
if (mapFlag[srcI])
{
const pointField
pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
const point& pt = pts[0];
pointIndexHit hit = tgt_.cellTree().findNearest(pt, GREAT);
if (hit.hit())
{
srcSeedI = srcI;
tgtSeedI = hit.index();
return true;
}
else
{
FatalErrorIn
(
"bool Foam::mapNearestMethod::findInitialSeeds"
"("
"const labelList&, "
"const boolList&, "
"const label, "
"label&, "
"label&"
") const"
)
<< "Unable to find nearest target cell"
<< " for source cell " << srcI
<< " with centre "
<< srcCells[srcI].centre(srcPts, srcFaces)
<< abort(FatalError);
}
break;
}
}
if (debug)
{
Pout<< "could not find starting seed" << endl;
}
return false;
}
void Foam::mapNearestMethod::calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
)
{
List<DynamicList<label> > srcToTgt(src_.nCells());
List<DynamicList<label> > tgtToSrc(tgt_.nCells());
const scalarField& srcVc = src_.cellVolumes();
const scalarField& tgtVc = tgt_.cellVolumes();
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
do
{
// find nearest tgt cell
findNearestCell(src_, tgt_, srcCellI, tgtCellI);
// store src/tgt cell pair
srcToTgt[srcCellI].append(tgtCellI);
tgtToSrc[tgtCellI].append(srcCellI);
// mark source cell srcCellI and tgtCellI as matched
mapFlag[srcCellI] = false;
// accumulate intersection volume
V_ += srcVc[srcCellI];
// find new source cell
setNextNearestCells
(
startSeedI,
srcCellI,
tgtCellI,
mapFlag,
srcCellIDs
);
}
while (srcCellI >= 0);
// for the case of multiple source cells per target cell, select the
// nearest source cell only and discard the others
const vectorField& srcCc = src_.cellCentres();
const vectorField& tgtCc = tgt_.cellCentres();
forAll(tgtToSrc, targetCellI)
{
if (tgtToSrc[targetCellI].size() > 1)
{
const vector& tgtC = tgtCc[tgtCellI];
DynamicList<label>& srcCells = tgtToSrc[targetCellI];
label srcCellI = srcCells[0];
scalar d = magSqr(tgtC - srcCc[srcCellI]);
for (label i = 1; i < srcCells.size(); i++)
{
label srcI = srcCells[i];
scalar dNew = magSqr(tgtC - srcCc[srcI]);
if (dNew < d)
{
d = dNew;
srcCellI = srcI;
}
}
srcCells.clear();
srcCells.append(srcCellI);
}
}
// If there are more target cells than source cells, some target cells
// might not yet be mapped
forAll(tgtToSrc, tgtCellI)
{
if (tgtToSrc[tgtCellI].empty())
{
label srcCellI = findMappedSrcCell(tgtCellI, tgtToSrc);
findNearestCell(tgt_, src_, tgtCellI, srcCellI);
tgtToSrc[tgtCellI].append(srcCellI);
}
}
// transfer addressing into persistent storage
forAll(srcToTgtCellAddr, i)
{
scalar v = srcVc[i];
srcToTgtCellAddr[i].transfer(srcToTgt[i]);
srcToTgtCellWght[i] = scalarList(1, v);
}
forAll(tgtToSrcCellAddr, i)
{
scalar v = tgtVc[i];
tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
tgtToSrcCellWght[i] = scalarList(1, v);
}
}
void Foam::mapNearestMethod::findNearestCell
(
const polyMesh& mesh1,
const polyMesh& mesh2,
const label cell1,
label& cell2
) const
{
const vectorField& Cc1 = mesh1.cellCentres();
const vectorField& Cc2 = mesh2.cellCentres();
const vector& p1 = Cc1[cell1];
DynamicList<label> cells2(10);
cells2.append(cell2);
DynamicList<label> visitedCells(10);
scalar d = GREAT;
do
{
label c2 = cells2.remove();
visitedCells.append(c2);
scalar dTest = magSqr(Cc2[c2] - p1);
if (dTest < d)
{
cell2 = c2;
d = dTest;
appendNbrCells(cell2, mesh2, visitedCells, cells2);
}
} while (cells2.size() > 0);
}
void Foam::mapNearestMethod::setNextNearestCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
boolList& mapFlag,
const labelList& srcCellIDs
) const
{
const labelList& srcNbr = src_.cellCells()[srcCellI];
srcCellI = -1;
forAll(srcNbr, i)
{
label cellI = srcNbr[i];
if (mapFlag[cellI])
{
srcCellI = cellI;
startSeedI = cellI + 1;
return;
}
}
(void)findInitialSeeds
(
srcCellIDs,
mapFlag,
startSeedI,
srcCellI,
tgtCellI
);
}
Foam::label Foam::mapNearestMethod::findMappedSrcCell
(
const label tgtCellI,
const List<DynamicList<label> >& tgtToSrc
) const
{
DynamicList<label> testCells(10);
DynamicList<label> visitedCells(10);
testCells.append(tgtCellI);
do
{
// search target tgtCellI neighbours for match with source cell
label tgtI = testCells.remove();
if (findIndex(visitedCells, tgtI) == -1)
{
visitedCells.append(tgtI);
if (tgtToSrc[tgtI].size())
{
return tgtToSrc[tgtI][0];
}
else
{
const labelList& nbrCells = tgt_.cellCells()[tgtI];
forAll(nbrCells, i)
{
if (findIndex(visitedCells, nbrCells[i]) == -1)
{
testCells.append(nbrCells[i]);
}
}
}
}
} while (testCells.size());
// did not find any match - should not be possible to get here!
return -1;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::mapNearestMethod::mapNearestMethod
(
const polyMesh& src,
const polyMesh& tgt
)
:
meshToMeshMethod(src, tgt)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::mapNearestMethod::~mapNearestMethod()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::mapNearestMethod::calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToSrcAddr,
scalarListList& tgtToSrcWght
)
{
bool ok = initialise
(
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght
);
if (!ok)
{
return;
}
// (potentially) participating source mesh cells
const labelList srcCellIDs(maskCells());
// list to keep track of whether src cell can be mapped
boolList mapFlag(src_.nCells(), false);
UIndirectList<bool>(mapFlag, srcCellIDs) = true;
// find initial point in tgt mesh
label srcSeedI = -1;
label tgtSeedI = -1;
label startSeedI = 0;
bool startWalk =
findInitialSeeds
(
srcCellIDs,
mapFlag,
startSeedI,
srcSeedI,
tgtSeedI
);
if (startWalk)
{
calculateAddressing
(
srcToTgtAddr,
srcToTgtWght,
tgtToSrcAddr,
tgtToSrcWght,
srcSeedI,
tgtSeedI,
srcCellIDs,
mapFlag,
startSeedI
);
}
else
{
// if meshes are collocated, after inflating the source mesh bounding
// box tgt mesh cells may be transferred, but may still not overlap
// with the source mesh
return;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,152 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::mapNearestMethod
Description
Map nearest mesh-to-mesh interpolation class
Not volume conservative.
SourceFiles
mapNearestMethod.C
\*---------------------------------------------------------------------------*/
#ifndef mapNearestMethod_H
#define mapNearestMethod_H
#include "meshToMeshMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mapNearestMethod Declaration
\*---------------------------------------------------------------------------*/
class mapNearestMethod
:
public meshToMeshMethod
{
protected:
// Protected Member Functions
//- Find indices of overlapping cells in src and tgt meshes - returns
// true if found a matching pair
bool findInitialSeeds
(
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const;
//- Calculate the mesh-to-mesh addressing and weights
void calculateAddressing
(
labelListList& srcToTgtCellAddr,
scalarListList& srcToTgtCellWght,
labelListList& tgtToSrcCellAddr,
scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
);
//- Find the nearest cell on mesh2 for cell1 on mesh1
void findNearestCell
(
const polyMesh& mesh1,
const polyMesh& mesh2,
const label cell1,
label& cell2
) const;
//- Set the next cells for the marching front algorithm
void setNextNearestCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
boolList& mapFlag,
const labelList& srcCellIDs
) const;
//- Find a source cell mapped to target cell tgtCellI
label findMappedSrcCell
(
const label tgtCellI,
const List<DynamicList<label> >& tgtToSrc
) const;
//- Disallow default bitwise copy construct
mapNearestMethod(const mapNearestMethod&);
//- Disallow default bitwise assignment
void operator=(const mapNearestMethod&);
public:
//- Run-time type information
TypeName("mapNearest");
//- Construct from source and target meshes
mapNearestMethod(const polyMesh& src, const polyMesh& tgt);
//- Destructor
virtual ~mapNearestMethod();
// Member Functions
// Evaluate
//- Calculate addressing and weights
virtual void calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToTgtAddr,
scalarListList& tgtToTgtWght
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,274 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "meshToMeshMethod.H"
#include "tetOverlapVolume.H"
#include "OFstream.H"
#include "Time.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(meshToMeshMethod, 0);
defineRunTimeSelectionTable(meshToMeshMethod, components);
}
Foam::scalar Foam::meshToMeshMethod::tolerance_ = 1e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::labelList Foam::meshToMeshMethod::maskCells() const
{
boundBox intersectBb
(
max(src_.bounds().min(), tgt_.bounds().min()),
min(src_.bounds().max(), tgt_.bounds().max())
);
intersectBb.inflate(0.01);
const cellList& srcCells = src_.cells();
const faceList& srcFaces = src_.faces();
const pointField& srcPts = src_.points();
DynamicList<label> cells(src_.nCells());
forAll(srcCells, srcI)
{
boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
if (intersectBb.overlaps(cellBb))
{
cells.append(srcI);
}
}
if (debug)
{
Pout<< "participating source mesh cells: " << cells.size() << endl;
}
return cells;
}
bool Foam::meshToMeshMethod::intersect
(
const label srcCellI,
const label tgtCellI
) const
{
scalar threshold = tolerance_*src_.cellVolumes()[srcCellI];
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell
(
pointField
(
tgt_.points(),
tgt_.cellPoints()[tgtCellI]
)
);
return overlapEngine.cellCellOverlapMinDecomp
(
src_,
srcCellI,
tgt_,
tgtCellI,
bbTgtCell,
threshold
);
}
Foam::scalar Foam::meshToMeshMethod::interVol
(
const label srcCellI,
const label tgtCellI
) const
{
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell
(
pointField
(
tgt_.points(),
tgt_.cellPoints()[tgtCellI]
)
);
scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
(
src_,
srcCellI,
tgt_,
tgtCellI,
bbTgtCell
);
return vol;
}
void Foam::meshToMeshMethod::appendNbrCells
(
const label cellI,
const polyMesh& mesh,
const DynamicList<label>& visitedCells,
DynamicList<label>& nbrCellIDs
) const
{
const labelList& nbrCells = mesh.cellCells()[cellI];
// filter out cells already visited from cell neighbours
forAll(nbrCells, i)
{
label nbrCellI = nbrCells[i];
if
(
(findIndex(visitedCells, nbrCellI) == -1)
&& (findIndex(nbrCellIDs, nbrCellI) == -1)
)
{
nbrCellIDs.append(nbrCellI);
}
}
}
bool Foam::meshToMeshMethod::initialise
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToSrcAddr,
scalarListList& tgtToSrcWght
) const
{
srcToTgtAddr.setSize(src_.nCells());
srcToTgtWght.setSize(src_.nCells());
tgtToSrcAddr.setSize(tgt_.nCells());
tgtToSrcWght.setSize(tgt_.nCells());
if (!src_.nCells())
{
return false;
}
else if (!tgt_.nCells())
{
if (debug)
{
Pout<< "mesh interpolation: hhave " << src_.nCells() << " source "
<< " cells but no target cells" << endl;
}
return false;
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::meshToMeshMethod::meshToMeshMethod
(
const polyMesh& src,
const polyMesh& tgt
)
:
src_(src),
tgt_(tgt),
V_(0.0)
{
if (!src_.nCells() || !tgt_.nCells())
{
if (debug)
{
Pout<< "mesh interpolation: cells not on processor: Source cells = "
<< src_.nCells() << ", target cells = " << tgt_.nCells()
<< endl;
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::meshToMeshMethod::~meshToMeshMethod()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::meshToMeshMethod::writeConnectivity
(
const polyMesh& mesh1,
const polyMesh& mesh2,
const labelListList& mesh1ToMesh2Addr
) const
{
Pout<< "Source size = " << mesh1.nCells() << endl;
Pout<< "Target size = " << mesh2.nCells() << endl;
word fName("addressing_" + mesh1.name() + "_to_" + mesh2.name());
if (Pstream::parRun())
{
fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
}
OFstream os(src_.time().path()/fName + ".obj");
label vertI = 0;
forAll(mesh1ToMesh2Addr, i)
{
const labelList& addr = mesh1ToMesh2Addr[i];
forAll(addr, j)
{
label cellI = addr[j];
const vector& c0 = mesh1.cellCentres()[i];
const cell& c = mesh2.cells()[cellI];
const pointField pts(c.points(mesh2.faces(), mesh2.points()));
forAll(pts, j)
{
const point& p = pts[j];
os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
vertI++;
os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
<< nl;
vertI++;
os << "l " << vertI - 1 << ' ' << vertI << nl;
}
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,188 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::meshToMeshMethod
Description
Base class for mesh-to-mesh calculation methods
SourceFiles
meshToMeshMethod.C
\*---------------------------------------------------------------------------*/
#ifndef meshToMeshMethod_H
#define meshToMeshMethod_H
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class meshToMeshMethod Declaration
\*---------------------------------------------------------------------------*/
class meshToMeshMethod
{
protected:
// Protected data
//- Reference to the source mesh
const polyMesh& src_;
//- Reference to the target mesh
const polyMesh& tgt_;
//- Cell total volume in overlap region [m3]
scalar V_;
//- Tolerance used in volume overlap calculations
static scalar tolerance_;
// Protected Member Functions
//- Return src cell IDs for the overlap region
labelList maskCells() const;
//- Return the true if cells intersect
bool intersect(const label srcCellI, const label tgtCellI) const;
//- Return the intersection volume between two cells
scalar interVol(const label srcCellI, const label tgtCellI) const;
//- Append target cell neihgbour cells to cellIDs list
void appendNbrCells
(
const label tgtCellI,
const polyMesh& mesh,
const DynamicList<label>& visitedTgtCells,
DynamicList<label>& nbrTgtCellIDs
) const;
bool initialise
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToTgtAddr,
scalarListList& tgtToTgtWght
) const;
//- Disallow default bitwise copy construct
meshToMeshMethod(const meshToMeshMethod&);
//- Disallow default bitwise assignment
void operator=(const meshToMeshMethod&);
public:
//- Run-time type information
TypeName("meshToMeshMethod");
//- Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
meshToMeshMethod,
components,
(
const polyMesh& src,
const polyMesh& tgt
),
(src, tgt)
);
//- Construct from source and target meshes
meshToMeshMethod(const polyMesh& src, const polyMesh& tgt);
//- Selector
static autoPtr<meshToMeshMethod> New
(
const word& methodName,
const polyMesh& src,
const polyMesh& tgt
);
//- Destructor
virtual ~meshToMeshMethod();
// Member Functions
// Evaluate
//- Calculate addressing and weights
virtual void calculate
(
labelListList& srcToTgtAddr,
scalarListList& srcToTgtWght,
labelListList& tgtToTgtAddr,
scalarListList& tgtToTgtWght
) = 0;
// Access
//- Return const access to the source mesh
inline const polyMesh& src() const;
//- Return const access to the target mesh
inline const polyMesh& tgt() const;
//- Return const access to the overlap volume
inline scalar V() const;
// Check
//- Write the connectivity (debugging)
void writeConnectivity
(
const polyMesh& mesh1,
const polyMesh& mesh2,
const labelListList& mesh1ToMesh2Addr
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "meshToMeshMethodI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
const Foam::polyMesh& Foam::meshToMeshMethod::src() const
{
return src_;
}
const Foam::polyMesh& Foam::meshToMeshMethod::tgt() const
{
return tgt_;
}
Foam::scalar Foam::meshToMeshMethod::V() const
{
return V_;
}
// ************************************************************************* //

View File

@ -0,0 +1,65 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "meshToMeshMethod.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::meshToMeshMethod> Foam::meshToMeshMethod::New
(
const word& methodName,
const polyMesh& src,
const polyMesh& tgt
)
{
if (debug)
{
Info<< "Selecting AMIMethod " << methodName << endl;
}
componentsConstructorTable::iterator cstrIter =
componentsConstructorTablePtr_->find(methodName);
if (cstrIter == componentsConstructorTablePtr_->end())
{
FatalErrorIn
(
"Foam::autoPtr<Foam::meshToMeshMethod> Foam::meshToMeshMethod::New"
"("
"const word&, "
"const polyMesh&, "
"const polyMesh&"
")"
) << "Unknown meshToMesh type "
<< methodName << nl << nl
<< "Valid meshToMesh types are:" << nl
<< componentsConstructorTablePtr_->sortedToc() << exit(FatalError);
}
return autoPtr<meshToMeshMethod>(cstrIter()(src, tgt));
}
// ************************************************************************* //

View File

@ -24,14 +24,9 @@ License
\*---------------------------------------------------------------------------*/
#include "meshToMeshNew.H"
#include "OFstream.H"
#include "Time.H"
#include "globalIndex.H"
#include "mergePoints.H"
#include "treeBoundBox.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "ListOps.H"
#include "meshToMeshMethod.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -55,55 +50,9 @@ namespace Foam
meshToMeshNew::interpolationMethodNames_;
}
Foam::scalar Foam::meshToMeshNew::tolerance_ = 1e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::meshToMeshNew::writeConnectivity
(
const polyMesh& src,
const polyMesh& tgt,
const labelListList& srcToTargetAddr
) const
{
Pout<< "Source size = " << src.nCells() << endl;
Pout<< "Target size = " << tgt.nCells() << endl;
word fName("addressing_" + src.name() + "_to_" + tgt.name());
if (Pstream::parRun())
{
fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
}
OFstream os(src.time().path()/fName + ".obj");
label vertI = 0;
forAll(srcToTargetAddr, i)
{
const labelList& tgtAddress = srcToTargetAddr[i];
forAll(tgtAddress, j)
{
label tgtI = tgtAddress[j];
const vector& c0 = src.cellCentres()[i];
const cell& c = tgt.cells()[tgtI];
const pointField pts(c.points(tgt.faces(), tgt.points()));
forAll(pts, j)
{
const point& p = pts[j];
os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
vertI++;
os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
<< nl;
vertI++;
os << "l " << vertI - 1 << ' ' << vertI << nl;
}
}
}
}
Foam::labelList Foam::meshToMeshNew::maskCells
(
const polyMesh& src,
@ -141,82 +90,6 @@ Foam::labelList Foam::meshToMeshNew::maskCells
}
bool Foam::meshToMeshNew::findInitialSeeds
(
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const
{
const cellList& srcCells = src.cells();
const faceList& srcFaces = src.faces();
const pointField& srcPts = src.points();
for (label i = startSeedI; i < srcCellIDs.size(); i++)
{
label srcI = srcCellIDs[i];
if (mapFlag[srcI])
{
const pointField
pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
forAll(pts, ptI)
{
const point& pt = pts[ptI];
label tgtI = tgt.cellTree().findInside(pt);
if (tgtI != -1 && intersect(src, tgt, srcI, tgtI))
{
srcSeedI = srcI;
tgtSeedI = tgtI;
return true;
}
}
}
}
if (debug)
{
Pout<< "could not find starting seed" << endl;
}
return false;
}
void Foam::meshToMeshNew::appendNbrCells
(
const label cellI,
const polyMesh& mesh,
const DynamicList<label>& visitedCells,
DynamicList<label>& nbrCellIDs
) const
{
const labelList& nbrCells = mesh.cellCells()[cellI];
// filter out cells already visited from cell neighbours
forAll(nbrCells, i)
{
label nbrCellI = nbrCells[i];
if
(
(findIndex(visitedCells, nbrCellI) == -1)
&& (findIndex(nbrCellIDs, nbrCellI) == -1)
)
{
nbrCellIDs.append(nbrCellI);
}
}
}
void Foam::meshToMeshNew::normaliseWeights
(
const word& descriptor,
@ -260,124 +133,29 @@ void Foam::meshToMeshNew::calcAddressing
const polyMesh& tgt
)
{
srcToTgtCellAddr_.setSize(src.nCells());
srcToTgtCellWght_.setSize(src.nCells());
tgtToSrcCellAddr_.setSize(tgt.nCells());
tgtToSrcCellWght_.setSize(tgt.nCells());
if (!src.nCells() || !tgt.nCells())
{
if (debug)
{
Pout<< "mesh interpolation: cells not on processor: Source cells = "
<< src.nCells() << ", target cells = " << tgt.nCells()
<< endl;
}
}
if (!src.nCells())
{
return;
}
else if (!tgt.nCells())
{
if (debug)
{
Pout<< "mesh interpolation: hhave " << src.nCells() << " source "
<< " cells but no target cells" << endl;
}
return;
}
// (potentially) participating source mesh cells
const labelList srcCellIDs = maskCells(src, tgt);
// list to keep track of whether src cell can be mapped
boolList mapFlag(src.nCells(), false);
UIndirectList<bool>(mapFlag, srcCellIDs) = true;
// find initial point in tgt mesh
label srcSeedI = -1;
label tgtSeedI = -1;
label startSeedI = 0;
bool startWalk =
findInitialSeeds
autoPtr<meshToMeshMethod> methodPtr
(
meshToMeshMethod::New
(
interpolationMethodNames_[method_],
src,
tgt,
srcCellIDs,
mapFlag,
startSeedI,
srcSeedI,
tgtSeedI
);
if (!startWalk)
{
// if meshes are collocated, after inflating the source mesh bounding
// box tgt mesh cells may be transferred, but may still not overlap
// with the source mesh
return;
}
switch (method_)
{
case imDirect:
{
calcDirect(src, tgt, srcSeedI, tgtSeedI);
break;
}
case imMapNearest:
{
calcMapNearest
(
src,
tgt,
srcSeedI,
tgtSeedI,
srcCellIDs,
mapFlag,
startSeedI
);
break;
}
case imCellVolumeWeight:
{
calcCellVolumeWeight
(
src,
tgt,
srcSeedI,
tgtSeedI,
srcCellIDs,
mapFlag,
startSeedI
);
break;
}
default:
{
FatalErrorIn
(
"void Foam::meshToMeshNew::calcAddressing"
"("
"const polyMesh&, "
"const polyMesh&"
")"
tgt
)
<< "Unknown interpolation method"
<< abort(FatalError);
}
}
);
methodPtr->calculate
(
srcToTgtCellAddr_,
srcToTgtCellWght_,
tgtToSrcCellAddr_,
tgtToSrcCellWght_
);
V_ = methodPtr->V();
if (debug)
{
writeConnectivity(src, tgt, srcToTgtCellAddr_);
methodPtr->writeConnectivity(src, tgt, srcToTgtCellAddr_);
}
}
@ -577,11 +355,59 @@ void Foam::meshToMeshNew::calculate()
}
Foam::AMIPatchToPatchInterpolation::interpolationMethod
Foam::meshToMeshNew::interpolationMethodAMI
(
const interpolationMethod method
) const
{
switch (method_)
{
case imDirect:
{
return AMIPatchToPatchInterpolation::imDirect;
break;
}
case imMapNearest:
{
return AMIPatchToPatchInterpolation::imMapNearest;
break;
}
case imCellVolumeWeight:
{
return AMIPatchToPatchInterpolation::imFaceAreaWeight;
break;
}
default:
{
FatalErrorIn
(
"Foam::AMIPatchToPatchInterpolation::interpolationMethod"
"Foam::meshToMeshNew::interpolationMethodAMI"
"("
"const interpolationMethod method"
") const"
)
<< "Unhandled enumeration " << method_
<< abort(FatalError);
}
}
return AMIPatchToPatchInterpolation::imDirect;
}
const Foam::PtrList<Foam::AMIPatchToPatchInterpolation>&
Foam::meshToMeshNew::patchAMIs() const
{
if (patchAMIs_.empty())
{
const word amiMethod =
AMIPatchToPatchInterpolation::interpolationMethodToWord
(
interpolationMethodAMI(method_)
);
patchAMIs_.setSize(srcPatchID_.size());
forAll(srcPatchID_, i)
@ -593,7 +419,9 @@ Foam::meshToMeshNew::patchAMIs() const
const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI];
Info<< "Creating AMI between source patch " << srcPP.name()
<< " and target patch " << tgtPP.name() << endl;
<< " and target patch " << tgtPP.name()
<< " using " << amiMethod
<< endl;
Info<< incrIndent;
@ -605,6 +433,7 @@ Foam::meshToMeshNew::patchAMIs() const
srcPP,
tgtPP,
faceAreaIntersect::tmMesh,
interpolationMethodAMI(method_),
true // flip target patch since patch normals are aligned
)
);

View File

@ -27,20 +27,14 @@ Class
Description
Class to calculate the cell-addressing between two overlapping meshes
Three methods are currently available:
- direct : 1-to-1 mapping between meshes
- mapNearest : assign nearest cell values without interpolation
- cellVolumeWeight : volume consistent mapping
The \c direct and \c cellVolumeWeight options are volume conservative,
whereas mapNearest is non-conservative.
Mapping is performed using a run-time selectable interpolation mothod
SeeAlso
meshToMeshMethod
SourceFiles
calcDirect.C
calcMapNearest.C
calcCellVolumeWeight.C
meshToMeshNew.C
meshToMeshNewParallelOps.C
meshToMeshNewTemplates.C
\*---------------------------------------------------------------------------*/
@ -70,7 +64,7 @@ public:
// Public data types
//- Enumeration specifying required accuracy
//- Enumeration specifying interpolation method
enum interpolationMethod
{
imDirect,
@ -128,9 +122,6 @@ private:
//- Target map pointer - parallel running only
autoPtr<mapDistribute> tgtMapPtr_;
//- Tolerance used in volume overlap calculations
static scalar tolerance_;
// Private Member Functions
@ -138,39 +129,9 @@ private:
template<class Type>
void add(UList<Type>& fld, const label offset) const;
//- Write the connectivity (debugging)
void writeConnectivity
(
const polyMesh& src,
const polyMesh& tgt,
const labelListList& srcToTargetAddr
) const;
//- Return src cell IDs for the overlap region
labelList maskCells(const polyMesh& src, const polyMesh& tgt) const;
//- Find indices of overlapping cells in src and tgt meshes - returns
// true if found a matching pair
bool findInitialSeeds
(
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs,
const boolList& mapFlag,
const label startSeedI,
label& srcSeedI,
label& tgtSeedI
) const;
//- Append target cell neihgbour cells to cellIDs list
void appendNbrCells
(
const label tgtCellI,
const polyMesh& tgt,
const DynamicList<label>& visitedTgtCells,
DynamicList<label>& nbrTgtCellIDs
) const;
//- Normalise the interpolation weights
void normaliseWeights
(
@ -180,121 +141,6 @@ private:
scalarListList& wght
) const;
// Direct (one-to-one) mapping
//- Main driver routine for direct mapping
void calcDirect
(
const polyMesh& src,
const polyMesh& tgt,
const label srcSeedI,
const label tgtSeedI
);
//- Append to list of src mesh seed indices
void appendToDirectSeeds
(
const polyMesh& src,
const polyMesh& tgt,
boolList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
label& srcSeedI,
label& tgtSeedI
) const;
// Nearest (non-conformal) mapping
//- Main driver routine for nearest-mapping routine
void calcMapNearest
(
const polyMesh& src,
const polyMesh& tgt,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
);
//- Find target cell index of cell closest to source cell
void findNearestCell
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
label& tgtCellI
);
//- Set the next pair of cells
void setNextNearestCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
boolList& mapFlag,
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs
);
//- Find source cell for target cell
label findMappedSrcCell
(
const polyMesh& tgt,
const label tgtCellI,
const List<DynamicList<label> >& tgtToSrc
) const;
// Cell volume weighted (non-conformal) interpolation
//- Main driver routine for cell volume weighted interpolation
void calcCellVolumeWeight
(
const polyMesh& src,
const polyMesh& tgt,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
boolList& mapFlag,
label& startSeedI
);
//- Set the next cells in the advancing front algorithm
void setNextCells
(
label& startSeedI,
label& srcCellI,
label& tgtCellI,
const polyMesh& src,
const polyMesh& tgt,
const labelList& srcCellIDs,
const boolList& mapFlag,
const DynamicList<label>& visitedCells,
labelList& seedCells
) const;
//- Return the true if cells intersect
bool intersect
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
const label tgtCellI
) const;
//- Return the intersection volume between two cells
scalar interVol
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
const label tgtCellI
) const;
//- Calculate the addressing between overalping regions of src and tgt
// meshes
void calcAddressing(const polyMesh& src, const polyMesh& tgt);
@ -302,6 +148,13 @@ private:
//- Calculate - main driver function
void calculate();
//- Conversion between mesh and patch interpolation methods
AMIPatchToPatchInterpolation::interpolationMethod
interpolationMethodAMI
(
const interpolationMethod method
) const;
//- Return the list of AMIs between source and target patches
const PtrList<AMIPatchToPatchInterpolation>& patchAMIs() const;

View File

@ -168,18 +168,32 @@ void Foam::sampledSurfaces::sampleAndWrite
template<class GeoField>
void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& objects)
{
wordList names;
if (loadFromFiles_)
{
IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
wordList names(fieldObjects.names());
labelList fieldNames(findStrings(fieldSelection_, names));
forAll(fieldNames, fieldI)
names = fieldObjects.names();
}
else
{
const word& fieldName = names[fieldNames[fieldI]];
names = mesh_.thisDb().names<GeoField>();
}
labelList nameIDs(findStrings(fieldSelection_, names));
wordHashSet fieldNames(wordList(names, nameIDs));
forAllConstIter(wordHashSet, fieldNames, iter)
{
const word& fieldName = iter.key();
if ((Pstream::master()) && verbose_)
{
Pout<< "sampleAndWrite: " << fieldName << endl;
}
if (loadFromFiles_)
{
const GeoField fld
(
IOobject
@ -192,30 +206,10 @@ void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& objects)
mesh_
);
if ((Pstream::master()) && verbose_)
{
Pout<< "sampleAndWrite: " << fieldName << endl;
}
sampleAndWrite(fld);
}
}
else
{
const wordList fieldNames
(
mesh_.thisDb().names<GeoField>(fieldSelection_)
);
forAll(fieldNames, i)
{
const word& fieldName = fieldNames[i];
if ((Pstream::master()) && verbose_)
{
Pout<< "sampleAndWrite: " << fieldName << endl;
}
sampleAndWrite
(
mesh_.thisDb().lookupObject<GeoField>(fieldName)

View File

@ -48,7 +48,8 @@ namespace Foam
defineTemplateTypeNameAndDebugWithName \
( \
sChemistryl##Comp##SThermo, \
(#sChemistry"<"#Comp"," + SThermo::typeName() + ">").c_str(), \
(word(sChemistry::typeName_()) + "<"#Comp"," + SThermo::typeName() + \
+ ">").c_str(), \
0 \
); \
\

View File

@ -47,8 +47,8 @@ namespace Foam
defineTemplateTypeNameAndDebugWithName \
( \
SS##Schem##Comp##SThermo##GThermo, \
(#SS"<"#Schem"<"#Comp"," + SThermo::typeName() + "," \
+ GThermo::typeName() + ">>").c_str(), \
(#SS"<" + word(Schem::typeName_()) +"<"#Comp"," + SThermo::typeName() \
+ "," + GThermo::typeName() + ">>").c_str(), \
0 \
); \
\

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -175,6 +175,7 @@ class triSurface
void writeDXGeometry(const bool, Ostream&) const;
void writeDXTrailer(Ostream&) const;
// Static private functions
//- Convert faces to labelledTri. All get same region.
@ -203,6 +204,7 @@ class triSurface
//- read non-comment line
static string getLineNoComment(IFstream&);
protected:
// Protected Member Functions
@ -219,6 +221,7 @@ protected:
return static_cast<List<Face>&>(*this);
}
public:
// Public typedefs
@ -287,7 +290,6 @@ public:
triSurface(const triSurface&);
//- Destructor
~triSurface();
@ -324,6 +326,7 @@ public:
// If >2 neighbours: undetermined.
const labelList& edgeOwner() const;
// Edit
//- Move points
@ -381,6 +384,7 @@ public:
labelList& faceMap
) const;
// Write
//- Write to Ostream in simple FOAM format

View File

@ -47,7 +47,7 @@ void dynLagrangian::updateSubGridScaleFields
const tmp<volTensorField>& gradU
)
{
nuSgs_ = (flm_/fmm_)*delta()*sqrt(k(gradU));
nuSgs_ = (flm_/fmm_)*sqr(delta())*mag(dev(symm(gradU)));
nuSgs_.correctBoundaryConditions();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -153,7 +153,10 @@ public:
//- Return SGS kinetic energy
tmp<volScalarField> k(const tmp<volTensorField>& gradU) const
{
return 2.0*sqr(delta())*magSqr(dev(symm(gradU)));
return
pow(2.0*flm_/fmm_, 2.0/3.0)
* pow(ce_, -2.0/3.0)
* sqr(delta())*magSqr(dev(symm(gradU)));
}
//- Return SGS kinetic energy

View File

@ -18,7 +18,7 @@ FoamFile
chemistryType
{
chemistrySolver ode;
chemistryThermo pyrolysisChemistryModel;
chemistryThermo pyrolysis;
}
chemistry on;