diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C new file mode 100644 index 0000000000000000000000000000000000000000..03f6355870b2b89e37940ad64b89bb0f48d1e1d7 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C @@ -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() +{} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H new file mode 100644 index 0000000000000000000000000000000000000000..7b6482b7c6cfe1850461b39895f088000f4e3119 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H @@ -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 + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodI.H new file mode 100644 index 0000000000000000000000000000000000000000..6a86eccf9c6f383fbef90ccf0ca857297098d3f7 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodI.H @@ -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_; +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodNew.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodNew.C new file mode 100644 index 0000000000000000000000000000000000000000..f2ba968316df07673b6e8c745ab9dcd0bfa208dd --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodNew.C @@ -0,0 +1,79 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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> +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 +) +{ + 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 + ) + ); +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C new file mode 100644 index 0000000000000000000000000000000000000000..9a47f132878fc4f99a683f5da6ead4174e419011 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C @@ -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(srcAddr[i].size(), magSf); + } + forAll(tgtAddr, i) + { + scalar magSf = this->tgtMagSf_[i]; + tgtAddress[i].transfer(tgtAddr[i]); + tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf); + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H new file mode 100644 index 0000000000000000000000000000000000000000..e7842ff9f8957b2557bf9245a6166b608655f38d --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H @@ -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 + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C new file mode 100644 index 0000000000000000000000000000000000000000..0356b52c0e670a7c0260cafe6129b9cd30ccff15 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C @@ -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]); + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H new file mode 100644 index 0000000000000000000000000000000000000000..0d76fe4d1bb785ec388abce9febef701d7dd6e16 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H @@ -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("faceAreaWeight"); + + + // 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 + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C new file mode 100644 index 0000000000000000000000000000000000000000..39e059273d15bb15ae18479f665fac92ca494c7c --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C @@ -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(srcAddr[i].size(), magSf); + } + forAll(tgtAddr, i) + { + scalar magSf = this->tgtMagSf_[i]; + tgtAddress[i].transfer(tgtAddr[i]); + tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf); + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H new file mode 100644 index 0000000000000000000000000000000000000000..0daf26662d27a98473c4b2e233ed405f65fa2aba --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H @@ -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 + +// ************************************************************************* //