diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index e281fc989f..d7864f7605 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -326,9 +326,30 @@ void Foam::AMIInterpolation::agglomerate // So now we have agglomeration of the target side in // allRestrict: - // 0..size-1 : local agglomeration (= targetRestrictAddressing) + // 0..size-1 : local agglomeration (= targetRestrictAddressing + // (but potentially permutated)) // size.. : agglomeration data from other processors + + // The trickiness in this algorithm is finding out the compaction + // of the remote data (i.e. allocation of the coarse 'slots'). We could + // either send across the slot compaction maps or just make sure + // that we allocate the slots in exactly the same order on both sending + // and receiving side (e.g. if the submap is set up to send 4 items, + // the constructMap is also set up to receive 4 items. + + + // Short note about the various types of indices: + // - face indices : indices into the geometry. + // - coarse face indices : how the faces get agglomerated + // - transferred data : how mapDistribute sends/receives data, + // - slots : indices into data after distribution (e.g. stencil, + // srcAddress/tgtAddress). Note: for fully local addressing + // the slots are equal to face indices. + // A mapDistribute has: + // - a subMap : these are face indices + // - a constructMap : these are from 'transferred-date' to slots + labelListList tgtSubMap(Pstream::nProcs()); // Local subMap is just identity @@ -340,10 +361,14 @@ void Foam::AMIInterpolation::agglomerate { if (procI != Pstream::myProcNo()) { - // Combine entries that point to the same coarse element. All - // the elements refer to local data so index into - // targetRestrictAddressing or allRestrict (since the same - // for local data). + // Combine entries that point to the same coarse element. + // The important bit is to loop over the data (and hand out + // compact indices ) in 'transferred data' order. This + // guarantees that we're doing exactly the + // same on sending and receiving side - e.g. the fourth element + // in the subMap is the fourth element received in the + // constructMap + const labelList& elems = map.subMap()[procI]; const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()]; @@ -383,8 +408,13 @@ void Foam::AMIInterpolation::agglomerate labelList tgtCompactMap(map.constructSize()); { - const labelList& elemsMap = - map.constructMap()[Pstream::myProcNo()]; + // Note that in special cases (e.g. 'appending' two AMIs) the + // local size after distributing can be longer than the number + // of faces. I.e. it duplicates elements. + // Since we don't know this size instead we loop over all + // reachable elements (using the local constructMap) + + const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()]; forAll(elemsMap, i) { label fineElem = elemsMap[i]; @@ -503,7 +533,7 @@ void Foam::AMIInterpolation::agglomerate forAll(fineSrcAddress, faceI) { // All the elements contributing to faceI. Are slots in - // mapDistribute'd data. + // target data. const labelList& elems = fineSrcAddress[faceI]; const scalarList& weights = fineSrcWeights[faceI]; const scalar fineArea = fineSrcMagSf[faceI]; diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C index 158fbce113..67bd613e1d 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C @@ -25,6 +25,8 @@ License #include "cyclicPeriodicAMIPolyPatch.H" #include "addToRunTimeSelectionTable.H" +#include "OBJstream.H" +#include "Time.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -263,6 +265,25 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI nbrPoints ); + + autoPtr thisStr; + autoPtr nbrStr; + if (debug) + { + const fileName dir(boundaryMesh().mesh().time().timePath()); + + mkDir(dir); + + thisStr.reset(new OBJstream(dir/name()+"_this.obj")); + nbrStr.reset(new OBJstream(dir/name()+"_nbr.obj")); + Info<< name() << " : dumping local duplicated geometry to " + << thisStr().name() << " and " << nbrStr().name() << endl; + + thisStr().write(thisPatch0, thisPatch0.points()); + nbrStr().write(nbrPatch0, nbrPatch0.points()); + } + + // Construct a new AMI interpolation between the initial patch locations AMIPtr_.reset ( @@ -287,12 +308,16 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI scalar srcSum(gAverage(AMIPtr_->srcWeightsSum())); scalar tgtSum(gAverage(AMIPtr_->tgtWeightsSum())); - // Direction of geometry replication + // Direction (or rather side of AMI : this or nbr patch) of + // geometry replication bool direction = nTransforms_ >= 0; // Increase in the source weight sum for the last iteration in the // opposite direction. If the current increase is less than this, the - // direction is reversed. + // direction (= side of AMI to transform) is reversed. + // We switch the side to replicate instead of reversing the transform + // since at the coupledPolyPatch level there is no + // 'reverseTransformPosition' functionality. scalar srcSumDiff = 0; // Loop, replicating the geometry @@ -311,6 +336,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI thisPatch.movePoints(thisPoints); + if (thisStr.valid()) + { + thisStr().write(thisPatch, thisPatch.points()); + } + AMIPtr_->append(thisPatch, nbrPatch0); } else @@ -319,6 +349,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI nbrPatch.movePoints(nbrPoints); + if (nbrStr.valid()) + { + nbrStr().write(nbrPatch, nbrPatch.points()); + } + AMIPtr_->append(thisPatch0, nbrPatch); } @@ -346,7 +381,7 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI // Check that the match is complete if (iter == maxIter_) { - FatalErrorIn + SeriousErrorIn ( "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI" "(" @@ -357,7 +392,12 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI << " do not couple to within a tolerance of " << matchTolerance() << " when transformed according to the periodic patch " - << periodicPatch.name() << "." << exit(FatalError); + << periodicPatch.name() << "." << nl + << "This is only acceptable during post-processing" + << "; not during running. Improve your mesh or increase" + << " the 'matchTolerance' setting in the patch specification." + << endl; + // << exit(FatalError); } // Check that both patches have replicated an integer number of times @@ -367,7 +407,7 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI || mag(tgtSum - floor(tgtSum + 0.5)) > tgtSum*matchTolerance() ) { - FatalErrorIn + SeriousErrorIn ( "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI" "(" @@ -377,7 +417,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI << "Patches " << name() << " and " << neighbPatch().name() << " do not overlap an integer number of times when transformed" << " according to the periodic patch " - << periodicPatch.name() << "." << exit(FatalError); + << periodicPatch.name() << "." << nl + << "This is only acceptable during post-processing" + << "; not during running. Improve your mesh or increase" + << " the 'matchTolerance' setting in the patch specification." + << endl; //<< exit(FatalError); } // Normalise the weights