ENH: Added initial support for particle tracking across AMI patches

src/lagrangian/basic/Cloud/Cloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,16 +44,22 @@ void Foam::Cloud<ParticleType>::checkPatches() const
     {
         if (isA<cyclicAMIPolyPatch>(pbm[patchI]))
         {
-            ok = false;
-            break;
+            const cyclicAMIPolyPatch& cami =
+                refCast<const cyclicAMIPolyPatch>(pbm[patchI]);
+
+            if (cami.owner())
+            {
+                ok = ok && (cami.AMI().singlePatchProc() != -1);
+            }
         }
     }
 
     if (!ok)
     {
-        WarningIn("void Foam::Cloud<ParticleType>::initCloud(const bool)")
-            << "Particle tracking across AMI patches is not currently "
-            << "supported" << endl;
+        FatalErrorIn("void Foam::Cloud<ParticleType>::initCloud(const bool)")
+            << "Particle tracking across AMI patches is only currently "
+            << "supported for cases where the AMI patches reside on a "
+            << "single processor" << abort(FatalError);
     }
 }
 

src/lagrangian/basic/particle/particle.H

@@ -257,6 +257,15 @@ protected:
         template<class TrackData>
         void hitCyclicPatch(const cyclicPolyPatch&, TrackData& td);
 
+        //- Overridable function to handle the particle hitting a cyclicAMIPatch
+        template<class TrackData>
+        void hitCyclicAMIPatch
+        (
+            const cyclicAMIPolyPatch&,
+            TrackData& td,
+            const vector& direction
+        );
+
         //- Overridable function to handle the particle hitting a
         //  processorPatch
         template<class TrackData>

src/lagrangian/basic/particle/particleTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "IOPosition.H"
 
 #include "cyclicPolyPatch.H"
+#include "cyclicAMIPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "wallPolyPatch.H"
@@ -598,6 +599,15 @@ Foam::scalar Foam::particle::trackToFace
                     static_cast<const cyclicPolyPatch&>(patch), td
                 );
             }
+            else if (isA<cyclicAMIPolyPatch>(patch))
+            {
+                p.hitCyclicAMIPatch
+                (
+                    static_cast<const cyclicAMIPolyPatch&>(patch),
+                    td,
+                    endPosition - position_
+                );
+            }
             else if (isA<processorPolyPatch>(patch))
             {
                 p.hitProcessorPatch
@@ -1010,6 +1020,78 @@ void Foam::particle::hitCyclicPatch
 }
 
 
+template<class TrackData>
+void Foam::particle::hitCyclicAMIPatch
+(
+    const cyclicAMIPolyPatch& cpp,
+    TrackData& td,
+    const vector& direction
+)
+{
+    const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
+
+    // patch face index on sending side
+    label patchFaceI = faceI_ - cpp.start();
+
+    // patch face index on receiving side - also updates position
+    patchFaceI = cpp.pointFace(patchFaceI, direction, position_);
+
+    if (patchFaceI < 0)
+    {
+        FatalErrorIn
+        (
+            "template<class TrackData>"
+            "void Foam::particle::hitCyclicAMIPatch"
+            "("
+                "const cyclicAMIPolyPatch&, "
+                "TrackData&, "
+                "const vector&"
+            ")"
+        )
+            << "Particle lost across " << cyclicAMIPolyPatch::typeName
+            << " patches " << cpp.name() << " and " << receiveCpp.name()
+            << " at position " << position_ << abort(FatalError);
+    }
+
+    // convert face index into global numbering
+    faceI_ = patchFaceI + receiveCpp.start();
+
+    cellI_ = mesh_.faceOwner()[faceI_];
+
+    tetFaceI_ = faceI_;
+
+    // See note in correctAfterParallelTransfer for tetPtI_ addressing.
+    tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
+
+    // Now the particle is on the receiving side
+
+    // Have patch transform the position
+    receiveCpp.transformPosition(position_, patchFaceI);
+
+    // Transform the properties
+    if (!receiveCpp.parallel())
+    {
+        const tensor& T =
+        (
+            receiveCpp.forwardT().size() == 1
+          ? receiveCpp.forwardT()[0]
+          : receiveCpp.forwardT()[patchFaceI]
+        );
+        transformProperties(T);
+    }
+    else if (receiveCpp.separated())
+    {
+        const vector& s =
+        (
+            (receiveCpp.separation().size() == 1)
+          ? receiveCpp.separation()[0]
+          : receiveCpp.separation()[patchFaceI]
+        );
+        transformProperties(-s);
+    }
+}
+
+
 template<class TrackData>
 void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&)
 {}

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -1248,6 +1248,114 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
 }
 
 
+template<class SourcePatch, class TargetPatch>
+Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcPointFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const vector& n,
+    const label tgtFaceI,
+    point& tgtPoint
+)
+const
+{
+    const pointField& srcPoints = srcPatch.points();
+
+    // source face addresses that intersect target face tgtFaceI
+    const labelList& addr = tgtAddress_[tgtFaceI];
+
+    forAll(addr, i)
+    {
+        label srcFaceI = addr[i];
+        const face& f = srcPatch[tgtFaceI];
+
+        pointHit ray = f.ray(tgtPoint, n, srcPoints);
+
+        if (ray.hit())
+        {
+            tgtPoint = ray.rawPoint();
+
+            return srcFaceI;
+        }
+    }
+
+    // no hit registered - try with face normal instead of input normal
+    forAll(addr, i)
+    {
+        label srcFaceI = addr[i];
+        const face& f = srcPatch[tgtFaceI];
+
+        vector nFace(-srcPatch.faceNormals()[srcFaceI]);
+        nFace += tgtPatch.faceNormals()[tgtFaceI];
+
+        pointHit ray = f.ray(tgtPoint, nFace, srcPoints);
+
+        if (ray.hit())
+        {
+            tgtPoint = ray.rawPoint();
+
+            return srcFaceI;
+        }
+    }
+
+    return -1;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtPointFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const vector& n,
+    const label srcFaceI,
+    point& srcPoint
+)
+const
+{
+    const pointField& tgtPoints = tgtPatch.points();
+
+    // target face addresses that intersect source face srcFaceI
+    const labelList& addr = srcAddress_[srcFaceI];
+
+    forAll(addr, i)
+    {
+        label tgtFaceI = addr[i];
+        const face& f = tgtPatch[tgtFaceI];
+
+        pointHit ray = f.ray(srcPoint, n, tgtPoints);
+
+        if (ray.hit())
+        {
+            srcPoint = ray.rawPoint();
+
+            return tgtFaceI;
+        }
+    }
+
+    // no hit registered - try with face normal instead of input normal
+    forAll(addr, i)
+    {
+        label tgtFaceI = addr[i];
+        const face& f = tgtPatch[tgtFaceI];
+
+        vector nFace(-srcPatch.faceNormals()[srcFaceI]);
+        nFace += tgtPatch.faceNormals()[tgtFaceI];
+
+        pointHit ray = f.ray(srcPoint, n, tgtPoints);
+
+        if (ray.hit())
+        {
+            srcPoint = ray.rawPoint();
+
+            return tgtFaceI;
+        }
+    }
+
+    return -1;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::writeFaceConnectivity
 (

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -454,6 +454,31 @@ public:
             ) const;
 
 
+        // Point intersections
+
+            //- Return source patch face index of point on target patch face
+            label srcPointFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const vector& n,
+                const label tgtFaceI,
+                point& tgtPoint
+            )
+            const;
+
+            //- Return target patch face index of point on source patch face
+            label tgtPointFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const vector& n,
+                const label srcFaceI,
+                point& srcPoint
+            )
+            const;
+
+
         // Checks
 
             //- Write face connectivity as OBJ file

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -797,6 +797,38 @@ bool Foam::cyclicAMIPolyPatch::order
 }
 
 
+Foam::label Foam::cyclicAMIPolyPatch::pointFace
+(
+    const label faceI,
+    const vector& n,
+    point& p
+) const
+{
+    if (owner())
+    {
+        return AMI().tgtPointFace
+        (
+            *this,
+            neighbPatch(),
+            n,
+            faceI,
+            p
+        );
+    }
+    else
+    {
+        return neighbPatch().AMI().srcPointFace
+        (
+            neighbPatch(),
+            *this,
+            n,
+            faceI,
+            p
+        );
+    }
+}
+
+
 void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
 {
     coupledPolyPatch::write(os);

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -367,6 +367,15 @@ public:
             labelList& rotation
         ) const;
 
+        //- Return face index on neighbour patch which shares point p
+        //  following trajectory vector n
+        label pointFace
+        (
+            const label faceI,
+            const vector& n,
+            point& p
+        ) const;
+
         //- Write the polyPatch data as a dictionary
         virtual void write(Ostream&) const;
 };