diff --git a/src/AMIInterpolation/AMIInterpolation.C b/src/AMIInterpolation/AMIInterpolation.C new file mode 100644 index 0000000000..0e5c6348e5 --- /dev/null +++ b/src/AMIInterpolation/AMIInterpolation.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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 +void Foam::AMIInterpolation::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 +void Foam::AMIInterpolation::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::checkPatches" + "(" + "const primitivePatch&, " + "const primitivePatch&" + ")" + ) << "Source patch is larger than, or misaligned with the " + << "projection surface" << endl; + } + + if (!bbSurf.contains(bbTgt)) + { + WarningIn + ( + "AMIInterpolation::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::checkPatches" + "(" + "const primitivePatch&, " + "const primitivePatch&" + ")" + ) << "Source and target patch bounding boxes are not similar" + << endl; + } +} + + +template +void Foam::AMIInterpolation::projectPointsToSurface +( + const searchableSurface& surf, + pointField& pts +) const +{ + List 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 +Foam::label Foam::AMIInterpolation::findTargetFace +( + const label srcFaceI, + const primitivePatch& srcPatch, + const primitivePatch& tgtPatch +) const +{ + label targetFaceI = -1; + + Random rndGen(123456); + + treeBoundBox bb(tgtPatch.points()); + + typedef treeDataPrimitivePatch treeType; + + indexedOctree 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 +void Foam::AMIInterpolation::appendNbrFaces +( + const label faceI, + const primitivePatch& patch, + const DynamicList