ENH: Updated AMI functionality - re: projection surface

This commit is contained in:
andy
2011-09-01 15:45:20 +01:00
parent 59ca09f152
commit 3a2b4dca63
5 changed files with 151 additions and 145 deletions

View File

@ -499,11 +499,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
pointField& pts
) const
{
if (!projectPoints_)
{
return;
}
if (debug)
{
Info<< "AMI: projecting points to surface" << endl;
@ -803,7 +798,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const searchableSurface& surf,
label srcFaceI,
label tgtFaceI
)
@ -824,6 +818,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
if (debug)
{
Info<< "AMI: calcAddressing" << endl;
writePatch(srcPatch, "VTK", "source");
writePatch(tgtPatch, "VTK", "target");
}
// temporary storage for addressing and weights
@ -833,57 +829,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
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("amiTgtPoints.obj");
forAll(tgtPoints, i)
{
meshTools::writeOBJ(os, tgtPoints[i]);
}
}
// map source and target patches onto projection surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
projectPointsToSurface(surf, srcPoints);
projectPointsToSurface(surf, tgtPoints);
// find initial face match using brute force/octree search
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if ((srcFaceI == -1) || (tgtFaceI == -1))
@ -891,9 +836,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
srcFaceI = 0;
tgtFaceI = 0;
bool foundFace = false;
forAll(srcPatch0, faceI)
forAll(srcPatch, faceI)
{
tgtFaceI = findTargetFace(srcFaceI, srcPatch0, tgtPatch0);
tgtFaceI = findTargetFace(faceI, srcPatch, tgtPatch);
if (tgtFaceI >= 0)
{
srcFaceI = faceI;
@ -911,7 +856,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
"("
"const primitivePatch&, "
"const primitivePatch&, "
"const searchableSurface&, "
"label, "
"label"
")"
@ -927,7 +871,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcPatch0.size();
label nFacesRemaining = srcPatch.size();
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
@ -953,7 +897,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
// append initial target face and neighbours
nbrFaces.append(tgtFaceI);
appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
bool faceProcessed = false;
@ -962,7 +906,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
// process new target face
tgtFaceI = nbrFaces.remove();
visitedFaces.append(tgtFaceI);
scalar area = interArea(srcFaceI, tgtFaceI, srcPatch0, tgtPatch0);
scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
// store when intersection area > 0
if (area > 0)
@ -973,7 +917,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
tgtAddr[tgtFaceI].append(srcFaceI);
tgtWght[tgtFaceI].append(area);
appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
faceProcessed = true;
}
@ -996,8 +940,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
(
srcFaceI,
tgtFaceI,
srcPatch0,
tgtPatch0,
srcPatch,
tgtPatch,
mapFlag,
seedFaces,
visitedFaces
@ -1084,9 +1028,8 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const searchableSurface& surf,
const faceAreaIntersect::triangulationMode& triMode,
const bool projectPoints
const autoPtr<searchableSurface>& surfPtr,
const faceAreaIntersect::triangulationMode& triMode
)
:
srcAddress_(),
@ -1095,24 +1038,74 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
tgtWeights_(),
startSeedI_(0),
triMode_(triMode),
projectPoints_(projectPoints),
srcMapPtr_(NULL),
tgtMapPtr_(NULL)
{
label srcSize = returnReduce(srcPatch.size(), sumOp<label>());
label tgtSize = returnReduce(tgtPatch.size(), sumOp<label>());
Info<< "AMI: Creating interpolation addressing and weights" << nl
<< " Source faces = " << srcSize << nl
<< " Target faces = " << tgtSize << endl;
Info<< "AMI: Creating addressing and weights between "
<< srcSize << " source faces and " << tgtSize << " target faces"
<< endl;
if (debug)
if (surfPtr.valid())
{
writePatch(srcPatch, "VTK", "source");
writePatch(tgtPatch, "VTK", "target");
}
// create new patches for source and target
pointField srcPoints = srcPatch.points();
primitivePatch srcPatch0
(
SubList<face>
(
srcPatch,
srcPatch.size(),
0
),
srcPoints
);
update(srcPatch, tgtPatch, surf);
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("amiTgtPoints.obj");
forAll(tgtPoints, i)
{
meshTools::writeOBJ(os, tgtPoints[i]);
}
}
// map source and target patches onto projection surface
projectPointsToSurface(surfPtr(), srcPoints);
projectPointsToSurface(surfPtr(), tgtPoints);
// calculate AMI interpolation
update(srcPatch0, tgtPatch0);
}
else
{
update(srcPatch, tgtPatch);
}
}
@ -1129,8 +1122,7 @@ template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const searchableSurface& surf
const primitivePatch& tgtPatch
)
{
static label patchI = 0;
@ -1180,7 +1172,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
// calculate AMI interpolation
calcAddressing(srcPatch, newTgtPatch, surf);
calcAddressing(srcPatch, newTgtPatch);
// Now
// ~~~
@ -1253,14 +1245,14 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
{
writeWeights(srcWeights_, srcPatch, "VTK", "source");
writeWeights(tgtWeights_, tgtPatch, "VTK", "target");
writeFaceConnectivity(srcPatch, newTgtPatch);
writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
}
}
else
{
checkPatches(srcPatch, tgtPatch);
calcAddressing(srcPatch, tgtPatch, surf);
calcAddressing(srcPatch, tgtPatch);
if (debug)
{
@ -1435,7 +1427,8 @@ template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::writeFaceConnectivity
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch
const primitivePatch& tgtPatch,
const labelListList& srcAddress
)
const
{
@ -1443,9 +1436,9 @@ const
label ptI = 1;
forAll(srcAddress_, i)
forAll(srcAddress, i)
{
const labelList& addr = srcAddress_[i];
const labelList& addr = srcAddress[i];
const point& srcPt = srcPatch.faceCentres()[i];
forAll(addr, j)
{

View File

@ -129,11 +129,10 @@ class AMIInterpolation
//- Face triangulation mode
const faceAreaIntersect::triangulationMode triMode_;
//- Flag to indicate whether or not to project points
const bool projectPoints_;
//- Source map pointer - parallel running only
autoPtr<mapDistribute> srcMapPtr_;
//- Target map pointer - parallel running only
autoPtr<mapDistribute> tgtMapPtr_;
@ -260,7 +259,6 @@ class AMIInterpolation
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const searchableSurface& surf,
label srcFaceI = -1,
label tgtFaceI = -1
);
@ -289,9 +287,8 @@ public:
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const searchableSurface& surf,
const faceAreaIntersect::triangulationMode& triMode,
const bool projectPoints
const autoPtr<searchableSurface>& surf,
const faceAreaIntersect::triangulationMode& triMode
);
@ -327,8 +324,7 @@ public:
void update
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const searchableSurface& surf
const primitivePatch& tgtPatch
);
@ -363,7 +359,8 @@ public:
void writeFaceConnectivity
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch
const primitivePatch& tgtPatch,
const labelListList& srcAddress
) const;
//- Write weights to VTK file

View File

@ -221,37 +221,12 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
}
const Foam::searchableSurface& Foam::cyclicAMIPolyPatch::surf()
void Foam::cyclicAMIPolyPatch::resetAMI()
{
word surfType(surfDict_.lookup("type"));
word surfName(surfDict_.lookupOrDefault("name", name()));
const polyMesh& mesh = boundaryMesh().mesh();
surfPtr_ =
searchableSurface::New
(
surfType,
IOobject
(
surfName,
mesh.time().constant(),
"triSurface",
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
surfDict_
);
return surfPtr_();
}
const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
{
if (!AMIPtr_.valid() && owner())
if (owner())
{
AMIPtr_.clear();
const polyPatch& nbr = nbrPatch();
pointField nbrPoints = nbrPatch().localPoints();
@ -282,8 +257,6 @@ const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
meshTools::writeOBJ(osO, this->localFaces(), localPoints());
}
bool projectPoints(readBool(surfDict_.lookup("projectPoints")));
// Construct/apply AMI interpolation to determine addressing and weights
AMIPtr_.reset
(
@ -291,14 +264,11 @@ const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
(
*this,
nbrPatch0,
surf(),
faceAreaIntersect::tmMesh,
projectPoints
surfPtr(),
faceAreaIntersect::tmMesh
)
);
}
return AMIPtr_();
}
@ -345,11 +315,7 @@ void Foam::cyclicAMIPolyPatch::movePoints
calcTransforms();
if (owner())
{
AMIPtr_.clear();
AMI();
}
resetAMI();
}
@ -611,6 +577,59 @@ bool Foam::cyclicAMIPolyPatch::owner() const
}
const Foam::autoPtr<Foam::searchableSurface>&
Foam::cyclicAMIPolyPatch::surfPtr()
{
word surfType(surfDict_.lookup("type"));
if (!surfPtr_.valid() && owner() && surfType != "none")
{
word surfName(surfDict_.lookupOrDefault("name", name()));
const polyMesh& mesh = boundaryMesh().mesh();
surfPtr_ =
searchableSurface::New
(
surfType,
IOobject
(
surfName,
mesh.time().constant(),
"triSurface",
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
surfDict_
);
}
return surfPtr_;
}
const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
{
if (!owner())
{
FatalErrorIn
(
"const AMIPatchToPatchInterpolation& cyclicAMIPolyPatch::AMI()"
)
<< "AMI interpolator only available to owner patch"
<< abort(FatalError);
}
if (!AMIPtr_.valid())
{
resetAMI();
}
return AMIPtr_();
}
void Foam::cyclicAMIPolyPatch::transformPosition(pointField& l) const
{
if (!parallel())
@ -704,11 +723,7 @@ void Foam::cyclicAMIPolyPatch::calcGeometry
nbrAreas
);
if (owner())
{
AMIPtr_.clear();
AMI();
}
resetAMI();
}

View File

@ -109,6 +109,9 @@ private:
const vectorField& half1Areas
);
//- Reset the AMI interpolator
void resetAMI();
protected:
@ -264,7 +267,7 @@ public:
inline const cyclicAMIPolyPatch& nbrPatch() const;
//- Return a reference to the projection surface
const searchableSurface& surf();
const autoPtr<searchableSurface>& surfPtr();
//- Return a reference to the AMI interpolator
const AMIPatchToPatchInterpolation& AMI();