mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: Initial commit of arbitrary mesh interface (AMI) functionality
This commit is contained in:
andy
793
src/AMIInterpolation/AMIInterpolation.C
Normal file
793
src/AMIInterpolation/AMIInterpolation.C
Normal file
@ -0,0 +1,793 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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 "AMIInterpolation.H"
|
||||
#include "faceAreaIntersect.H"
|
||||
#include "Random.H"
|
||||
#include "treeDataPrimitivePatch.H"
|
||||
#include "indexedOctree.H"
|
||||
#include "primitivePatch.H"
|
||||
#include "meshTools.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<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);
|
||||
|
||||
Info<< "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::AMIInterpolation<SourcePatch, TargetPatch>::checkPatches
|
||||
(
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
)
|
||||
{
|
||||
const scalar maxBoundsError = 0.05;
|
||||
|
||||
// sanity checks
|
||||
boundBox bbSrc(srcPatch.points());
|
||||
boundBox bbTgt(tgtPatch.points());
|
||||
boundBox bbSurf(srcPatch.points());
|
||||
|
||||
|
||||
// projection surface bounds - check against bounds of source and target
|
||||
|
||||
bbSurf.inflate(maxBoundsError);
|
||||
|
||||
if (!bbSurf.contains(bbSrc))
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"AMIInterpolation<SourcePatch, TargetPatch>::checkPatches"
|
||||
"("
|
||||
"const primitivePatch&, "
|
||||
"const primitivePatch&"
|
||||
")"
|
||||
) << "Source patch is larger than, or misaligned with the "
|
||||
<< "projection surface" << endl;
|
||||
}
|
||||
|
||||
if (!bbSurf.contains(bbTgt))
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"AMIInterpolation<SourcePatch, TargetPatch>::checkPatches"
|
||||
"("
|
||||
"const primitivePatch&, "
|
||||
"const primitivePatch&"
|
||||
")"
|
||||
) << "Target patch is larger than, or misaligned with the "
|
||||
<< "projection surface" << endl;
|
||||
}
|
||||
|
||||
|
||||
// check bounds of source and target
|
||||
|
||||
bbTgt.inflate(maxBoundsError);
|
||||
|
||||
if (!bbTgt.contains(bbSrc))
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"AMIInterpolation<SourcePatch, TargetPatch>::checkPatches"
|
||||
"("
|
||||
"const primitivePatch&, "
|
||||
"const primitivePatch&"
|
||||
")"
|
||||
) << "Source and target patch bounding boxes are not similar"
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
|
||||
(
|
||||
const searchableSurface& surf,
|
||||
pointField& pts
|
||||
) const
|
||||
{
|
||||
List<pointIndexHit> nearInfo;
|
||||
|
||||
surf.findNearest(pts, scalarField(pts.size(), GREAT), nearInfo);
|
||||
|
||||
label nMiss = 0;
|
||||
forAll(nearInfo, i)
|
||||
{
|
||||
const pointIndexHit& pi = nearInfo[i];
|
||||
|
||||
if (pi.hit())
|
||||
{
|
||||
pts[i] = pi.hitPoint();
|
||||
}
|
||||
else
|
||||
{
|
||||
pts[i] = pts[i];
|
||||
nMiss++;
|
||||
}
|
||||
}
|
||||
|
||||
if (nMiss > 0)
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void Foam::projectPointsToSurface"
|
||||
"("
|
||||
"const searchableSurface&, "
|
||||
"pointField&"
|
||||
") const"
|
||||
)
|
||||
<< "Error projecting points to surface: "
|
||||
<< nMiss << " faces could not be determined"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::findTargetFace
|
||||
(
|
||||
const label srcFaceI,
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
) const
|
||||
{
|
||||
label targetFaceI = -1;
|
||||
|
||||
Random rndGen(123456);
|
||||
|
||||
treeBoundBox bb(tgtPatch.points());
|
||||
|
||||
typedef treeDataPrimitivePatch<face, SubList, const pointField&> treeType;
|
||||
|
||||
indexedOctree<treeType> tree
|
||||
(
|
||||
treeType(false, tgtPatch),
|
||||
bb.extend(rndGen, 1E-4), // overall search domain
|
||||
8, // maxLevel
|
||||
10, // leafsize
|
||||
3.0 // duplicity
|
||||
);
|
||||
|
||||
const pointField& srcPts = srcPatch.points();
|
||||
const face& srcFace = srcPatch[srcFaceI];
|
||||
const point& srcPt = srcFace.centre(srcPts);
|
||||
const scalar srcFaceArea = srcFace.mag(srcPts);
|
||||
|
||||
// pointIndexHit sample = tree.findNearest(srcPt, sqr(0.1*bb.mag()));
|
||||
pointIndexHit sample = tree.findNearest(srcPt, 10.0*srcFaceArea);
|
||||
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Info<< "Source point = " << srcPt << ", Sample point = "
|
||||
<< sample.hitPoint() << ", Sample index = " << sample.index()
|
||||
<< endl;
|
||||
}
|
||||
|
||||
if (sample.hit())
|
||||
{
|
||||
targetFaceI = sample.index();
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"Foam::label Foam::cyclicAMIPolyPatch::findTargetFace"
|
||||
"("
|
||||
"const label, "
|
||||
"const primitivePatch&"
|
||||
"const primitivePatch&"
|
||||
") const"
|
||||
) << "Unable to find target face for source face" << srcFaceI
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
return targetFaceI;
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
|
||||
(
|
||||
const label faceI,
|
||||
const primitivePatch& patch,
|
||||
const DynamicList<label>& visitedFaces,
|
||||
DynamicList<label>& faceIDs
|
||||
) const
|
||||
{
|
||||
// const labelList& nbrFaces = patch.pointFaces()[faceI];
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
|
||||
(
|
||||
label& srcFaceI,
|
||||
label& tgtFaceI,
|
||||
const primitivePatch& srcPatch0,
|
||||
const primitivePatch& tgtPatch0,
|
||||
const boolList& mapFlag,
|
||||
labelList& seedFaces,
|
||||
const DynamicList<label>& visitedFaces
|
||||
) const
|
||||
{
|
||||
// const labelList& srcNbrFaces = srcPatch0.pointFaces()[srcFaceI];
|
||||
const labelList& srcNbrFaces = srcPatch0.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, srcPatch0, tgtPatch0);
|
||||
|
||||
if (area > 0)
|
||||
{
|
||||
// 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
|
||||
forAll(mapFlag, faceI)
|
||||
{
|
||||
if (mapFlag[faceI] && seedFaces[faceI] != -1)
|
||||
{
|
||||
srcFaceI = faceI;
|
||||
tgtFaceI = seedFaces[faceI];
|
||||
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// perform new search to find match
|
||||
if (debug)
|
||||
{
|
||||
Info<< "Advancing front stalled: searching for new "
|
||||
<< "target face" << endl;
|
||||
}
|
||||
|
||||
forAll(mapFlag, faceI)
|
||||
{
|
||||
if (mapFlag[faceI])
|
||||
{
|
||||
srcFaceI = faceI;
|
||||
tgtFaceI = findTargetFace(srcFaceI, srcPatch0, tgtPatch0);
|
||||
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
FatalErrorIn
|
||||
(
|
||||
"void Foam::cyclicAMIPolyPatch::setNextFaces"
|
||||
"("
|
||||
"label&, "
|
||||
"label&, "
|
||||
"const primitivePatch&, "
|
||||
"const primitivePatch&, "
|
||||
"const boolList&, "
|
||||
"const labelList&, "
|
||||
"const DynamicList<label>&"
|
||||
") const"
|
||||
) << "Unable to set source and target faces" << abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
|
||||
(
|
||||
const label srcFaceI,
|
||||
const label tgtFaceI,
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
) const
|
||||
{
|
||||
const pointField& srcPoints = srcPatch.points();
|
||||
const pointField& tgtPoints = tgtPatch.points();
|
||||
|
||||
const face& src = srcPatch[srcFaceI];
|
||||
const face& tgt = tgtPatch[tgtFaceI];
|
||||
|
||||
// quick reject if either face has zero area
|
||||
if ((src.mag(srcPoints) < ROOTVSMALL) || (tgt.mag(tgtPoints) < ROOTVSMALL))
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
// create intersection object
|
||||
faceAreaIntersect inter(srcPoints, tgtPoints);
|
||||
|
||||
// crude resultant norm
|
||||
const vector n = 0.5*(tgt.normal(tgtPoints) - src.normal(srcPoints));
|
||||
|
||||
scalar area = inter.calc(src, tgt, n);
|
||||
|
||||
if ((debug > 1) && (area > 0))
|
||||
{
|
||||
writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
|
||||
}
|
||||
|
||||
return area;
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
|
||||
(
|
||||
const word& patchType,
|
||||
const primitivePatch& patch,
|
||||
const pointField& points,
|
||||
const List<DynamicList<label> >& addr,
|
||||
List<DynamicList<scalar> >& wght
|
||||
)
|
||||
{
|
||||
scalarList wghtSum(patch.size(), 0.0);
|
||||
|
||||
// normalise weights by face areas
|
||||
forAll(wght, faceI)
|
||||
{
|
||||
const DynamicList<label>& addressing = addr[faceI];
|
||||
forAll(addressing, addrI)
|
||||
{
|
||||
const scalar faceArea = patch[faceI].mag(points);
|
||||
wght[faceI][addrI] /= faceArea;
|
||||
wghtSum[faceI] += wght[faceI][addrI];
|
||||
}
|
||||
}
|
||||
|
||||
Info<< "Cumulative " << patchType << "source weights min/max = "
|
||||
<< min(wghtSum) << ", " << max(wghtSum) << endl;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
|
||||
(
|
||||
const SourcePatch& srcPatch,
|
||||
const TargetPatch& tgtPatch,
|
||||
const searchableSurface& surf
|
||||
)
|
||||
:
|
||||
srcAddress_(),
|
||||
srcWeights_(),
|
||||
tgtAddress_(),
|
||||
tgtWeights_()
|
||||
{
|
||||
checkPatches(srcPatch, tgtPatch);
|
||||
|
||||
calcAddressing(srcPatch, tgtPatch, surf);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::~AMIInterpolation()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
|
||||
(
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch,
|
||||
const searchableSurface& surf
|
||||
)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "AMI: calcAddressing" << endl;
|
||||
}
|
||||
|
||||
// temporary storage for addressing and weights
|
||||
List<DynamicList<label> > srcAddr(srcPatch.size());
|
||||
List<DynamicList<scalar> > srcWght(srcPatch.size());
|
||||
List<DynamicList<label> > tgtAddr(tgtPatch.size());
|
||||
List<DynamicList<scalar> > tgtWght(tgtPatch.size());
|
||||
|
||||
|
||||
// create new patches for source and target
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
pointField srcPoints = srcPatch.points();
|
||||
primitivePatch srcPatch0
|
||||
(
|
||||
SubList<face>
|
||||
(
|
||||
srcPatch,
|
||||
srcPatch.size(),
|
||||
0
|
||||
),
|
||||
srcPoints
|
||||
);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
OFstream os("amiSrcPoints.obj");
|
||||
forAll(srcPoints, i)
|
||||
{
|
||||
meshTools::writeOBJ(os, srcPoints[i]);
|
||||
}
|
||||
}
|
||||
|
||||
pointField tgtPoints = tgtPatch.points();
|
||||
primitivePatch tgtPatch0
|
||||
(
|
||||
SubList<face>
|
||||
(
|
||||
tgtPatch,
|
||||
tgtPatch.size(),
|
||||
0
|
||||
),
|
||||
tgtPoints
|
||||
);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
OFstream os("osTgtPoints.obj");
|
||||
forAll(tgtPoints, i)
|
||||
{
|
||||
meshTools::writeOBJ(os, tgtPoints[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Info<< "AMI: projecting points to surface" << endl;
|
||||
}
|
||||
|
||||
// map source and target patches onto projection surface
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
projectPointsToSurface(surf, srcPoints);
|
||||
projectPointsToSurface(surf, tgtPoints);
|
||||
|
||||
|
||||
// find initial face match using brute force/octree search
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
label srcFaceI = 0;
|
||||
label tgtFaceI = findTargetFace(srcFaceI, srcPatch0, tgtPatch0);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Info<< "AMI: initial target face = " << tgtFaceI << endl;
|
||||
}
|
||||
|
||||
|
||||
// construct weights and addressing
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
label facesRemaining = srcPatch0.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(facesRemaining, -1);
|
||||
seedFaces[srcFaceI] = tgtFaceI;
|
||||
|
||||
// list to keep track of whether src face can be mapped
|
||||
boolList mapFlag(facesRemaining, true);
|
||||
|
||||
do
|
||||
{
|
||||
nbrFaces.clear();
|
||||
visitedFaces.clear();
|
||||
|
||||
// append initial target face and neighbours
|
||||
nbrFaces.append(tgtFaceI);
|
||||
appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
|
||||
|
||||
do
|
||||
{
|
||||
// process new target face
|
||||
tgtFaceI = nbrFaces.remove();
|
||||
visitedFaces.append(tgtFaceI);
|
||||
scalar area = interArea(srcFaceI, tgtFaceI, srcPatch0, tgtPatch0);
|
||||
|
||||
// store when intersection area > 0
|
||||
if (area > 0)
|
||||
{
|
||||
srcAddr[srcFaceI].append(tgtFaceI);
|
||||
srcWght[srcFaceI].append(area);
|
||||
|
||||
tgtAddr[tgtFaceI].append(srcFaceI);
|
||||
tgtWght[tgtFaceI].append(area);
|
||||
|
||||
appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
|
||||
}
|
||||
|
||||
} while (nbrFaces.size() > 0);
|
||||
|
||||
mapFlag[srcFaceI] = false;
|
||||
|
||||
facesRemaining--;
|
||||
|
||||
// choose new src face from current src face neighbour
|
||||
if (facesRemaining > 0)
|
||||
{
|
||||
setNextFaces
|
||||
(
|
||||
srcFaceI,
|
||||
tgtFaceI,
|
||||
srcPatch0,
|
||||
tgtPatch0,
|
||||
mapFlag,
|
||||
seedFaces,
|
||||
visitedFaces
|
||||
);
|
||||
}
|
||||
} while (facesRemaining > 0);
|
||||
|
||||
|
||||
// weights normalisation
|
||||
normaliseWeights("source", srcPatch0, srcPoints, srcAddr, srcWght);
|
||||
normaliseWeights("target", tgtPatch0, tgtPoints, tgtAddr, tgtWght);
|
||||
|
||||
// transfer data to persistent storage
|
||||
srcAddress_.setSize(srcAddr.size());
|
||||
srcWeights_.setSize(srcWght.size());
|
||||
forAll(srcAddr, i)
|
||||
{
|
||||
srcAddress_[i].transfer(srcAddr[i]);
|
||||
srcWeights_[i].transfer(srcWght[i]);
|
||||
}
|
||||
|
||||
tgtAddress_.setSize(tgtAddr.size());
|
||||
tgtWeights_.setSize(tgtWght.size());
|
||||
forAll(tgtAddr, i)
|
||||
{
|
||||
tgtAddress_[i].transfer(tgtAddr[i]);
|
||||
tgtWeights_[i].transfer(tgtWght[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::movePoints()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::Field<Type> >
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
|
||||
(
|
||||
const Field<Type>& fld
|
||||
) const
|
||||
{
|
||||
if (fld.size() != tgtAddress_.size())
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"AMIInterpolation::interpolateToSource(const Field<Type>)"
|
||||
) << "Supplied field size is not equal to target patch size. "
|
||||
<< "Target patch = " << tgtAddress_.size() << ", supplied field: "
|
||||
<< fld.size() << abort(FatalError);
|
||||
}
|
||||
|
||||
tmp<Field<Type> > tresult
|
||||
(
|
||||
new Field<Type>(srcAddress_.size(), pTraits<Type>::zero)
|
||||
);
|
||||
|
||||
Field<Type>& result = tresult();
|
||||
|
||||
forAll(result, faceI)
|
||||
{
|
||||
const labelList& faces = srcAddress_[faceI];
|
||||
const scalarList& weights = srcWeights_[faceI];
|
||||
|
||||
forAll(faces, i)
|
||||
{
|
||||
result[faceI] += fld[faces[i]]*weights[i];
|
||||
}
|
||||
}
|
||||
|
||||
return tresult;
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::Field<Type> >
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
|
||||
(
|
||||
const tmp<Field<Type> >& tFld
|
||||
) const
|
||||
{
|
||||
return interpolateToSource(tFld());
|
||||
}
|
||||
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::Field<Type> >
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
|
||||
(
|
||||
const Field<Type>& fld
|
||||
) const
|
||||
{
|
||||
if (fld.size() != srcAddress_.size())
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"AMIInterpolation::interpolateToSource(const Field<Type>)"
|
||||
) << "Supplied field size is not equal to source patch size. "
|
||||
<< "Source patch = " << srcAddress_.size() << ", supplied field: "
|
||||
<< fld.size() << abort(FatalError);
|
||||
}
|
||||
|
||||
tmp<Field<Type> > tresult
|
||||
(
|
||||
new Field<Type>(tgtAddress_.size(), pTraits<Type>::zero)
|
||||
);
|
||||
|
||||
Field<Type>& result = tresult();
|
||||
|
||||
forAll(result, faceI)
|
||||
{
|
||||
const labelList& faces = tgtAddress_[faceI];
|
||||
const scalarList& weights = tgtWeights_[faceI];
|
||||
|
||||
forAll(faces, i)
|
||||
{
|
||||
result[faceI] += fld[faces[i]]*weights[i];
|
||||
}
|
||||
}
|
||||
|
||||
return tresult;
|
||||
}
|
||||
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::Field<Type> >
|
||||
Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
|
||||
(
|
||||
const tmp<Field<Type> >& tFld
|
||||
) const
|
||||
{
|
||||
return interpolateToTarget(tFld());
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
258
src/AMIInterpolation/AMIInterpolation.H
Normal file
258
src/AMIInterpolation/AMIInterpolation.H
Normal file
@ -0,0 +1,258 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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::AMIInterpolation
|
||||
|
||||
Description
|
||||
Interpolation class dealing with transfer of data between two
|
||||
primitive patches with an arbitrary mesh interface (AMI)
|
||||
|
||||
SourceFiles
|
||||
AMIInterpolation.C
|
||||
AMIInterpolationName.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef AMIInterpolation_H
|
||||
#define AMIInterpolation_H
|
||||
|
||||
#include "className.H"
|
||||
#include "DynamicList.H"
|
||||
#include "searchableSurface.H"
|
||||
#include "boolList.H"
|
||||
#include "primitivePatch.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class AMIInterpolationName Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
TemplateName(AMIInterpolation);
|
||||
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class AMIInterpolation Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
template<class SourcePatch, class TargetPatch>
|
||||
class AMIInterpolation
|
||||
:
|
||||
public AMIInterpolationName
|
||||
{
|
||||
// Static data
|
||||
|
||||
// Source patch
|
||||
|
||||
//- Addressing
|
||||
labelListList srcAddress_;
|
||||
|
||||
//- Weighting factors
|
||||
scalarListList srcWeights_;
|
||||
|
||||
|
||||
// Target patch
|
||||
|
||||
//- Addressing
|
||||
labelListList tgtAddress_;
|
||||
|
||||
//- Weighting factors
|
||||
scalarListList tgtWeights_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Disallow default bitwise copy construct
|
||||
AMIInterpolation(const AMIInterpolation&);
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const AMIInterpolation&);
|
||||
|
||||
//- Write triangle intersection to OBJ file
|
||||
void writeIntersectionOBJ
|
||||
(
|
||||
const scalar area,
|
||||
const face& f1,
|
||||
const face& f2,
|
||||
const pointField& f1Points,
|
||||
const pointField& f2Points
|
||||
) const;
|
||||
|
||||
void checkPatches
|
||||
(
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
);
|
||||
|
||||
void projectPointsToSurface
|
||||
(
|
||||
const searchableSurface& surf,
|
||||
pointField& pts
|
||||
) const;
|
||||
|
||||
label findTargetFace
|
||||
(
|
||||
const label srcFaceI,
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
) const;
|
||||
|
||||
void appendNbrFaces
|
||||
(
|
||||
const label faceI,
|
||||
const primitivePatch& patch,
|
||||
const DynamicList<label>& visitedFaces,
|
||||
DynamicList<label>& faceIDs
|
||||
) const;
|
||||
|
||||
void setNextFaces
|
||||
(
|
||||
label& srcFaceI,
|
||||
label& tgtFaceI,
|
||||
const primitivePatch& srcPatch0,
|
||||
const primitivePatch& tgtPatch0,
|
||||
const boolList& mapFlag,
|
||||
labelList& seedFaces,
|
||||
const DynamicList<label>& visitedFaces
|
||||
) const;
|
||||
|
||||
scalar interArea
|
||||
(
|
||||
const label srcFaceI,
|
||||
const label tgtFaceI,
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch
|
||||
) const;
|
||||
|
||||
void normaliseWeights
|
||||
(
|
||||
const word& patchType,
|
||||
const primitivePatch& patch,
|
||||
const pointField& points,
|
||||
const List<DynamicList<label> >& addr,
|
||||
List<DynamicList<scalar> >& wght
|
||||
);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
AMIInterpolation
|
||||
(
|
||||
const SourcePatch& srcPatch,
|
||||
const TargetPatch& tgtPatch,
|
||||
const searchableSurface& surf
|
||||
);
|
||||
|
||||
|
||||
//- Destructor
|
||||
~AMIInterpolation();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Access
|
||||
|
||||
// Source patch
|
||||
|
||||
//- Return const access to source patch addressing
|
||||
inline const labelListList& srcAddress();
|
||||
|
||||
//- Return const access to source patch weights
|
||||
inline const scalarListList& srcWeights();
|
||||
|
||||
|
||||
// Target patch
|
||||
|
||||
//- Return const access to target patch addressing
|
||||
inline const labelListList& tgtAddress();
|
||||
|
||||
//- Return const access to target patch weights
|
||||
inline const scalarListList& tgtWeights();
|
||||
|
||||
|
||||
// Manipulation
|
||||
|
||||
//- Calculate addressing
|
||||
void calcAddressing
|
||||
(
|
||||
const primitivePatch& srcPatch,
|
||||
const primitivePatch& tgtPatch,
|
||||
const searchableSurface& surf
|
||||
);
|
||||
|
||||
//- Correct weighting factors for moving mesh
|
||||
bool movePoints();
|
||||
|
||||
|
||||
// Evaluation
|
||||
|
||||
//- Interpolate from target to source
|
||||
template<class Type>
|
||||
tmp<Field<Type> > interpolateToSource(const Field<Type>& fld) const;
|
||||
|
||||
//- Interpolate from target tmp field to source
|
||||
template<class Type>
|
||||
tmp<Field<Type> > interpolateToSource
|
||||
(
|
||||
const tmp<Field<Type> >& tFld
|
||||
) const;
|
||||
|
||||
//- Interpolate from source to target
|
||||
template<class Type>
|
||||
tmp<Field<Type> > interpolateToTarget(const Field<Type>& fld) const;
|
||||
|
||||
//- Interpolate from source tmp field to target
|
||||
template<class Type>
|
||||
tmp<Field<Type> > interpolateToTarget
|
||||
(
|
||||
const tmp<Field<Type> >& tFld
|
||||
) const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "AMIInterpolationI.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
# include "AMIInterpolation.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
32
src/AMIInterpolation/AMIInterpolationName.C
Normal file
32
src/AMIInterpolation/AMIInterpolationName.C
Normal file
@ -0,0 +1,32 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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 "AMIInterpolation.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
defineTypeNameAndDebug(Foam::AMIInterpolationName, 0);
|
||||
|
||||
// ************************************************************************* //
|
||||
50
src/AMIInterpolation/AMIPatchToPatchInterpolation.H
Normal file
50
src/AMIInterpolation/AMIPatchToPatchInterpolation.H
Normal file
@ -0,0 +1,50 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef AMIPatchToPatchInterpolation_H
|
||||
#define AMIPatchToPatchInterpolation_H
|
||||
|
||||
#include "AMIInterpolation.H"
|
||||
#include "PrimitivePatch.H"
|
||||
#include "face.H"
|
||||
#include "SubList.H"
|
||||
#include "pointField.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
typedef AMIInterpolation
|
||||
<
|
||||
PrimitivePatch<face, SubList, const pointField&>,
|
||||
PrimitivePatch<face, SubList, const pointField&>
|
||||
> AMIPatchToPatchInterpolation;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
4
src/AMIInterpolation/Make/files
Normal file
4
src/AMIInterpolation/Make/files
Normal file
@ -0,0 +1,4 @@
|
||||
AMIInterpolationName.C
|
||||
faceAreaIntersect/faceAreaIntersect.C
|
||||
|
||||
LIB = $(FOAM_USER_LIBBIN)/libAMIInterpolation
|
||||
6
src/AMIInterpolation/Make/options
Normal file
6
src/AMIInterpolation/Make/options
Normal file
@ -0,0 +1,6 @@
|
||||
EXE_INC = \
|
||||
-DFULLDEBUG -g -O0 \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude
|
||||
|
||||
LIB_LIBS = \
|
||||
-lmeshTools
|
||||
277
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
Normal file
277
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
Normal file
@ -0,0 +1,277 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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 "faceAreaIntersect.H"
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
void Foam::faceAreaIntersect::triSliceWithPlane
|
||||
(
|
||||
const triPoints& tri,
|
||||
const plane& p,
|
||||
FixedList<triPoints, 10>& tris,
|
||||
label& nTris
|
||||
)
|
||||
{
|
||||
// distance to cutting plane
|
||||
FixedList<scalar, 3> d;
|
||||
|
||||
// determine how many of the points are above the cutting plane
|
||||
label nCoPlanar = 0;
|
||||
label nPos = 0;
|
||||
label posI = -1;
|
||||
label negI = -1;
|
||||
bool coPlanar;
|
||||
forAll(tri, i)
|
||||
{
|
||||
d[i] = ((tri[i] - p.refPoint()) & p.normal());
|
||||
|
||||
if (mag(d[i]) < 1e-10)
|
||||
{
|
||||
coPlanar = true;
|
||||
nCoPlanar++;
|
||||
}
|
||||
else
|
||||
{
|
||||
coPlanar = false;
|
||||
}
|
||||
|
||||
if (!coPlanar)
|
||||
{
|
||||
if (d[i] > 0)
|
||||
{
|
||||
nPos++;
|
||||
posI = i;
|
||||
}
|
||||
else
|
||||
{
|
||||
negI = i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ((nPos == 3) || ((nPos == 2) && (nCoPlanar == 1)))
|
||||
{
|
||||
// all points above cutting plane - add triangle to list
|
||||
tris[nTris++] = tri;
|
||||
}
|
||||
else if ((nPos == 2) || ((nPos == 1) && (nCoPlanar == 1)))
|
||||
{
|
||||
// 2 points above plane, 1 below
|
||||
// resulting quad above plane split into 2 triangles
|
||||
|
||||
// point under the plane
|
||||
label i0 = negI;
|
||||
|
||||
// indices of remaining points
|
||||
label i1 = d.fcIndex(i0);
|
||||
label i2 = d.fcIndex(i1);
|
||||
|
||||
// determine the two intersection points
|
||||
point p01 = planeIntersection(d, tri, i0, i1);
|
||||
point p02 = planeIntersection(d, tri, i0, i2);
|
||||
|
||||
// forget triangle below plane
|
||||
// - decompose quad above plane into 2 triangles and add to list
|
||||
setTriPoints(tri[i1], tri[i2], p02, nTris, tris);
|
||||
setTriPoints(tri[i1], p02, p01, nTris, tris);
|
||||
}
|
||||
else if ((nPos == 1) && (nCoPlanar != 1))
|
||||
{
|
||||
// 1 point above plane, 2 below
|
||||
// resulting quad below plane split into 2 triangles
|
||||
|
||||
// point above the plane
|
||||
label i0 = posI;
|
||||
|
||||
// indices of remaining points
|
||||
label i1 = d.fcIndex(i0);
|
||||
label i2 = d.fcIndex(i1);
|
||||
|
||||
// determine the two intersection points
|
||||
point p01 = planeIntersection(d, tri, i0, i1);
|
||||
point p02 = planeIntersection(d, tri, i0, i2);
|
||||
|
||||
// forget quad below plane
|
||||
// - add triangle above plane to list
|
||||
setTriPoints(tri[i0], p01, p02, nTris, tris);
|
||||
}
|
||||
else
|
||||
{
|
||||
// all points below cutting plane - forget
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::faceAreaIntersect::triangleIntersect
|
||||
(
|
||||
const triPoints& src,
|
||||
const triPoints& tgt,
|
||||
const vector& n
|
||||
)
|
||||
{
|
||||
// Work storage
|
||||
FixedList<triPoints, 10> cutTris;
|
||||
label nCutTris = 0;
|
||||
|
||||
FixedList<triPoints, 10> tris;
|
||||
label nTris = 0;
|
||||
|
||||
// cut source triangle with all inwards pointing faces of target triangle
|
||||
// - triangles in cutTris are inside target triangle
|
||||
|
||||
// edge 0
|
||||
{
|
||||
// cut triangle src with plane and put resulting sub-triangles in
|
||||
// cutTris list
|
||||
|
||||
plane pl0(tgt[0], tgt[1], tgt[1] + n);
|
||||
triSliceWithPlane(src, pl0, cutTris, nCutTris);
|
||||
}
|
||||
|
||||
if (nCutTris == 0)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
// edge1
|
||||
{
|
||||
// cut cutTris with plane and put resulting sub-triangles in
|
||||
// tris list (re-use tris storage)
|
||||
|
||||
plane pl1(tgt[1], tgt[2], tgt[2] + n);
|
||||
|
||||
nTris = 0;
|
||||
|
||||
for (label i = 0; i < nCutTris; i++)
|
||||
{
|
||||
triSliceWithPlane(cutTris[i], pl1, tris, nTris);
|
||||
}
|
||||
|
||||
if (nTris == 0)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// edge2
|
||||
{
|
||||
// cut tris with plane and put resulting sub-triangles in
|
||||
// cutTris list (re-use cutTris storage)
|
||||
|
||||
plane pl2(tgt[2], tgt[0], tgt[0] + n);
|
||||
|
||||
nCutTris = 0;
|
||||
|
||||
for (label i = 0; i < nTris; i++)
|
||||
{
|
||||
triSliceWithPlane(tris[i], pl2, cutTris, nCutTris);
|
||||
}
|
||||
|
||||
if (nCutTris == 0)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
// calculate area of sub-triangles
|
||||
scalar area = 0.0;
|
||||
for (label i = 0; i < nCutTris; i++)
|
||||
{
|
||||
const triPoints& t = cutTris[i];
|
||||
area += mag(0.5*((t[1] - t[0])^(t[2] - t[0])));
|
||||
}
|
||||
|
||||
return area;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::faceAreaIntersect::faceAreaIntersect
|
||||
(
|
||||
const pointField& pointsA,
|
||||
const pointField& pointsB
|
||||
)
|
||||
:
|
||||
pointsA_(pointsA),
|
||||
pointsB_(pointsB)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::faceAreaIntersect::calc
|
||||
(
|
||||
const face& faceA,
|
||||
const face& faceB,
|
||||
const vector& n
|
||||
)
|
||||
{
|
||||
// split faces into triangles
|
||||
DynamicList<face> trisA;
|
||||
DynamicList<face> trisB;
|
||||
|
||||
// if (useTriangleFan)
|
||||
// {
|
||||
// triangleFan(faceA, trisA);
|
||||
// triangleFan(faceB, trisB);
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// faceA.triangles(pointsA_, trisA);
|
||||
// faceB.triangles(pointsB_, trisB);
|
||||
// }
|
||||
|
||||
faceA.triangles(pointsA_, trisA);
|
||||
faceB.triangles(pointsB_, trisB);
|
||||
|
||||
// intersect triangles
|
||||
scalar totalArea = 0.0;
|
||||
forAll(trisA, tA)
|
||||
{
|
||||
triPoints tpA = getTriPoints(pointsA_, trisA[tA], false);
|
||||
|
||||
if (triArea(tpA) > ROOTVSMALL)
|
||||
{
|
||||
forAll(trisB, tB)
|
||||
{
|
||||
triPoints tpB = getTriPoints(pointsB_, trisB[tB], true);
|
||||
|
||||
if (triArea(tpB) > ROOTVSMALL)
|
||||
{
|
||||
totalArea += triangleIntersect(tpA, tpB, n);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return totalArea;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
161
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
Normal file
161
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
Normal file
@ -0,0 +1,161 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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::faceAreaIntersect
|
||||
|
||||
Description
|
||||
Face intersection class
|
||||
- calculates intersection area by sub-dividing face into triangles
|
||||
and cutting
|
||||
|
||||
SourceFiles
|
||||
faceAreaIntersect.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef faceAreaIntersect_H
|
||||
#define faceAreaIntersect_H
|
||||
|
||||
#include "pointField.H"
|
||||
#include "FixedList.H"
|
||||
#include "plane.H"
|
||||
#include "face.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class faceAreaIntersect Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class faceAreaIntersect
|
||||
{
|
||||
typedef FixedList<point, 3> triPoints;
|
||||
|
||||
// Private data
|
||||
|
||||
//- Reference to the points of sideA
|
||||
const pointField& pointsA_;
|
||||
|
||||
//- Reference to the points of sideB
|
||||
const pointField& pointsB_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Get triPoints from face
|
||||
inline triPoints getTriPoints
|
||||
(
|
||||
const pointField& points,
|
||||
const face& f,
|
||||
const bool reverse
|
||||
) const;
|
||||
|
||||
//- Set triPoints into tri list
|
||||
inline void setTriPoints
|
||||
(
|
||||
const point& a,
|
||||
const point& b,
|
||||
const point& c,
|
||||
label& count,
|
||||
FixedList<triPoints, 10>& tris
|
||||
) const;
|
||||
|
||||
//- Decompose face into triangle fan
|
||||
inline void triangleFan
|
||||
(
|
||||
const face& f,
|
||||
DynamicList<face>& faces
|
||||
) const;
|
||||
|
||||
//- Return point of intersection between plane and triangle edge
|
||||
inline point planeIntersection
|
||||
(
|
||||
const FixedList<scalar, 3>& d,
|
||||
const triPoints& t,
|
||||
const label negI,
|
||||
const label posI
|
||||
) const;
|
||||
|
||||
//- Return triangle area
|
||||
inline scalar triArea(const triPoints& t) const;
|
||||
|
||||
|
||||
//- Slice triangle with plane and generate new cut sub-triangles
|
||||
void triSliceWithPlane
|
||||
(
|
||||
const triPoints& tri,
|
||||
const plane& p,
|
||||
FixedList<triPoints, 10>& tris,
|
||||
label& nTris
|
||||
);
|
||||
|
||||
//- Return area of intersection of triangles src and tgt
|
||||
scalar triangleIntersect
|
||||
(
|
||||
const triPoints& src,
|
||||
const triPoints& tgt,
|
||||
const vector& n
|
||||
);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
faceAreaIntersect
|
||||
(
|
||||
const pointField& pointsA,
|
||||
const pointField& pointsB
|
||||
);
|
||||
|
||||
|
||||
// Public Member Functions
|
||||
|
||||
//- Return area of intersection of faceA with faceB
|
||||
scalar calc
|
||||
(
|
||||
const face& faceA,
|
||||
const face& faceB,
|
||||
const vector& n
|
||||
);
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "faceAreaIntersectI.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
109
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
Normal file
109
src/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
Normal file
@ -0,0 +1,109 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 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/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
inline void Foam::faceAreaIntersect::setTriPoints
|
||||
(
|
||||
const point& a,
|
||||
const point& b,
|
||||
const point& c,
|
||||
label& count,
|
||||
FixedList<triPoints, 10>& tris
|
||||
) const
|
||||
{
|
||||
triPoints& tp = tris[count++];
|
||||
tp[0] = a;
|
||||
tp[1] = b;
|
||||
tp[2] = c;
|
||||
}
|
||||
|
||||
|
||||
inline Foam::faceAreaIntersect::triPoints Foam::faceAreaIntersect::getTriPoints
|
||||
(
|
||||
const pointField& points,
|
||||
const face& f,
|
||||
const bool reverse
|
||||
) const
|
||||
{
|
||||
triPoints result;
|
||||
|
||||
if (reverse)
|
||||
{
|
||||
result[2] = points[f[0]];
|
||||
result[1] = points[f[1]];
|
||||
result[0] = points[f[2]];
|
||||
}
|
||||
else
|
||||
{
|
||||
result[0] = points[f[0]];
|
||||
result[1] = points[f[1]];
|
||||
result[2] = points[f[2]];
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
inline void Foam::faceAreaIntersect::triangleFan
|
||||
(
|
||||
const face& f,
|
||||
DynamicList<face>& faces
|
||||
) const
|
||||
{
|
||||
if (f.size() > 2)
|
||||
{
|
||||
const label v0 = 0;
|
||||
|
||||
labelList indices(3);
|
||||
|
||||
for (label i = 1; i < f.size() - 1; i++)
|
||||
{
|
||||
indices[0] = f[v0];
|
||||
indices[1] = f[i];
|
||||
indices[2] = f[i + 1];
|
||||
faces.append(face(indices));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline Foam::point Foam::faceAreaIntersect::planeIntersection
|
||||
(
|
||||
const FixedList<scalar, 3>& d,
|
||||
const triPoints& t,
|
||||
const label negI,
|
||||
const label posI
|
||||
) const
|
||||
{
|
||||
return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
|
||||
{
|
||||
return mag(0.5*((t[1] - t[0])^(t[2] - t[0])));
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
|
||||
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -163,6 +163,21 @@ Foam::tmp<Foam::pointField> Foam::boundBox::points() const
|
||||
}
|
||||
|
||||
|
||||
void Foam::boundBox::inflate(const scalar s)
|
||||
{
|
||||
boundBox bb(*this);
|
||||
|
||||
vector dir = bb.span();
|
||||
scalar magDir = Foam::mag(dir);
|
||||
dir /= magDir;
|
||||
|
||||
scalar ext = s*magDir;
|
||||
|
||||
min_ -= dir*ext;
|
||||
max_ += dir*ext;
|
||||
}
|
||||
|
||||
|
||||
bool Foam::boundBox::contains(const UList<point>& points) const
|
||||
{
|
||||
if (points.empty())
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
|
||||
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -160,6 +160,13 @@ public:
|
||||
//- Return corner points in an order corresponding to a 'hex' cell
|
||||
tmp<pointField> points() const;
|
||||
|
||||
|
||||
// Manipulate
|
||||
|
||||
//- Inflate box by factor*mag(span) in all dimensions
|
||||
void inflate(const scalar s);
|
||||
|
||||
|
||||
// Query
|
||||
|
||||
//- Overlaps/touches boundingBox?
|
||||
|
||||
Reference in New Issue
Block a user