Merge branch 'HrvojeJasak'

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C

@@ -69,7 +69,7 @@ void Foam::cyclicGgiPolyPatch::checkDefinition() const
     // patches must specify the same information If not, well,
     // we stop the simulation and ask for a fix.
 
-    if (shadowIndex() < 0)
+    if (!active())
     {
         // No need to check anything, the shadow is not initialized properly.
         // This will happen with blockMesh when defining cyclicGGI patches.
@@ -243,7 +243,7 @@ const Foam::cyclicGgiPolyPatch& Foam::cyclicGgiPolyPatch::cyclicShadow() const
 
 void Foam::cyclicGgiPolyPatch::calcTransforms()
 {
-    if (debug)
+    if (active() && debug)
     {
         // Check definition of the cyclic pair
         checkDefinition();

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C

@@ -238,6 +238,15 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const
 
 void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
 {
+    if (reconFaceCellCentresPtr_)
+    {
+        FatalErrorIn
+        (
+            "void ggiPolyPatch::calcReconFaceCellCentres() const"
+        )   << "Reconstructed cell centres already calculated"
+            << abort(FatalError);
+    }
+
     // Create neighbouring face centres using interpolation
     if (master())
     {
@@ -277,19 +286,26 @@ void Foam::ggiPolyPatch::calcLocalParallel() const
     localParallelPtr_ = new bool(false);
     bool& emptyOrComplete = *localParallelPtr_;
 
-    // If running in serial, all GGIs are local parallel
+    // If running in serial, all GGIs are expanded to zone size.
+    // This happens on decomposition and reconstruction where
+    // size and shadow size may be zero, but zone size may not
     // HJ, 1/Jun/2011
     if (!Pstream::parRun())
     {
-        return;
+        emptyOrComplete = false;
     }
+    else
+    {
+        // Calculate localisation on master and shadow
+        emptyOrComplete =
+            (
+                zone().size() == size()
+             && shadow().zone().size() == shadow().size()
+            )
+         || (size() == 0 && shadow().size() == 0);
 
-    // Calculate localisation on master and shadow
+        reduce(emptyOrComplete, andOp<bool>());
-    emptyOrComplete =
+    }
-        (zone().size() == size() && shadow().zone().size() == shadow().size())
-     || (size() == 0 && shadow().size() == 0);
-
-    reduce(emptyOrComplete, andOp<bool>());
 
     if (debug && Pstream::parRun())
     {
@@ -607,9 +623,6 @@ Foam::label Foam::ggiPolyPatch::shadowIndex() const
         }
     }
 
-    // Force local parallel
-    localParallel();
-
     return shadowIndex_;
 }
 
@@ -720,12 +733,18 @@ void Foam::ggiPolyPatch::initAddressing()
 {
     if (active())
     {
-        // Force zone addressing first
+        // Calculate transforms for correct GGI cut
+        calcTransforms();
+
+        // Force zone addressing and remote zone addressing
+        // (uses GGI interpolator)
         zoneAddressing();
         remoteZoneAddressing();
 
+        // Force local parallel
         if (Pstream::parRun() && !localParallel())
         {
+            // Calculate send addressing
             sendAddr();
         }
     }
@@ -746,8 +765,6 @@ void Foam::ggiPolyPatch::initGeometry()
     // patch comms.  HJ, 11/Jul/2011
     if (active())
     {
-        calcTransforms();
-
         // Note: Only master calculates recon; slave uses master interpolation
         if (master())
         {

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H

@@ -108,9 +108,6 @@ class ggiPolyPatch
 
     // Private member functions
 
-        //- Is the GGI active& (zone and shadow present)
-        bool active() const;
-
         //- Calculate patch-to-zone addressing
         virtual void calcZoneAddressing() const;
 
@@ -159,6 +156,9 @@ protected:
 
     // Protected Member functions
 
+        //- Is the GGI active? (zone and shadow present)
+        bool active() const;
+
         //- Initialise the calculation of the patch addressing
         virtual void initAddressing();
 
@@ -285,15 +285,15 @@ public:
             return shadowName_;
         }
 
-        //- Return shadow patch index
-        label shadowIndex() const;
-
         //- Return name of interpolation face zone
         const word& zoneName() const
         {
             return zoneName_;
         }
 
+        //- Return shadow patch index
+        label shadowIndex() const;
+
         //- Return zone patch index
         label zoneIndex() const;
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C

@@ -313,6 +313,12 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::interpolate
     {
         // No expansion or filtering needed.  HJ, 4/Jun/2011
 
+        if (empty())
+        {
+            // Patch empty, no interpolation
+            return tmp<Field<Type> >(new Field<Type>());
+        }
+
         // Interpolate field
         if (master())
         {

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.C

@@ -49,29 +49,54 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::overlapGgiPolyPatch::active() const
+{
+    polyPatchID shadow(shadowName_, boundaryMesh());
+    faceZoneID zone(zoneName_, boundaryMesh().mesh().faceZones());
+
+    return shadow.active() && zone.active();
+}
+
+
 void Foam::overlapGgiPolyPatch::calcLocalParallel() const
 {
     // Calculate patch-to-zone addressing
     if (localParallelPtr_)
     {
-        FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const")
+        FatalErrorIn("void overlapGgiPolyPatch::calcLocalParallel() const")
             << "Local parallel switch already calculated"
             << abort(FatalError);
     }
 
+    // If running in serial, all GGIs are local parallel
+    // HJ, 1/Jun/2011
     localParallelPtr_ = new bool(false);
     bool& emptyOrComplete = *localParallelPtr_;
 
-    // Calculate localisation on master and shadow
+    // If running in serial, all GGIs are expanded to zone size.
-    emptyOrComplete =
+    // This happens on decomposition and reconstruction where
-        (zone().size() == size() && shadow().zone().size() == shadow().size())
+    // size and shadow size may be zero, but zone size may not
-     || (size() == 0 && shadow().size() == 0);
+    // HJ, 1/Jun/2011
+    if (!Pstream::parRun())
+    {
+        emptyOrComplete = false;
+    }
+    else
+    {
+        // Calculate localisation on master and shadow
+        emptyOrComplete =
+            (
+                zone().size() == size()
+             && shadow().zone().size() == shadow().size()
+            )
+         || (size() == 0 && shadow().size() == 0);
 
-    reduce(emptyOrComplete, andOp<bool>());
+        reduce(emptyOrComplete, andOp<bool>());
+    }
 
     if (debug && Pstream::parRun())
     {
-        Info<< "GGI patch Master: " << name()
+        Info<< "Overlap GGI patch Master: " << name()
             << " Slave: " << shadowName() << " is ";
 
         if (emptyOrComplete)
@@ -101,7 +126,7 @@ Foam::overlapGgiPolyPatch::overlapGgiPolyPatch
     zoneName_(word::null),
     shadowIndex_(-1),
     zoneIndex_(-1),
-    rotationAxis_(vector(0.0, 0.0, 1.0)),
+    rotationAxis_(vector(0, 0, 1)),
     nCopies_(0),
     expandedMasterPtr_(NULL),
     expandedSlavePtr_(NULL),
@@ -257,9 +282,6 @@ Foam::label Foam::overlapGgiPolyPatch::shadowIndex() const
         }
     }
 
-    // Force local parallel
-    localParallel();
-
     return shadowIndex_;
 }
 
@@ -290,11 +312,13 @@ Foam::overlapGgiPolyPatch::shadow() const
     return refCast<const overlapGgiPolyPatch>(boundaryMesh()[shadowIndex()]);
 }
 
+
 const Foam::faceZone& Foam::overlapGgiPolyPatch::zone() const
 {
     return boundaryMesh().mesh().faceZones()[zoneIndex()];
 }
 
+
 bool Foam::overlapGgiPolyPatch::localParallel() const
 {
     // Calculate patch-to-zone addressing
@@ -306,6 +330,7 @@ bool Foam::overlapGgiPolyPatch::localParallel() const
     return *localParallelPtr_;
 }
 
+
 const Foam::label& Foam::overlapGgiPolyPatch::nCopies() const
 {
     // Read the number of copies to be made from the dictionary for the
@@ -313,6 +338,7 @@ const Foam::label& Foam::overlapGgiPolyPatch::nCopies() const
     return nCopies_;
 }
 
+
 bool Foam::overlapGgiPolyPatch::master() const
 {
     // The first overlapggi interface is master,second one is slave

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H

@@ -147,6 +147,9 @@ protected:
 
     // Protected Member functions
 
+        //- Is the GGI active? (zone and shadow present)
+        bool active() const;
+
         //- Initialise the calculation of the patch addressing
         virtual void initAddressing();
 
@@ -301,7 +304,7 @@ public:
         //- Return wedge angle
         scalar angle() const
         {
-            return 360.0/nCopies();
+            return 360.0/scalar(nCopies());
         }
 
         //- Return number of slave copies

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatchGeometry.C

@@ -123,23 +123,24 @@ Foam::overlapGgiPolyPatch::calcExpandedGeometry(label ncp, label index) const
     return new standAlonePatch(expandedFaces, expandedPoints);
 }
 
+
 const Foam::standAlonePatch& Foam::overlapGgiPolyPatch::expandedMaster() const
 {
     if (!expandedMasterPtr_)
     {
-        expandedMasterPtr_ =
+        expandedMasterPtr_ = calcExpandedGeometry(nCopies(), zoneIndex());
-            calcExpandedGeometry( nCopies(), zoneIndex() );
     }
 
     return *expandedMasterPtr_;
 }
 
+
 const Foam::standAlonePatch& Foam::overlapGgiPolyPatch::expandedSlave() const
 {
     if (!expandedSlavePtr_)
     {
         expandedSlavePtr_ =
-            calcExpandedGeometry( shadow().nCopies(), shadow().zoneIndex() );
+            calcExpandedGeometry(shadow().nCopies(), shadow().zoneIndex());
     }
 
     return *expandedSlavePtr_;
@@ -168,7 +169,11 @@ void Foam::overlapGgiPolyPatch::calcPatchToPatch() const
                 reverseT(),
                 separation(),
                 0,             // master overlap tolerance
-                0              // slave overlap tolerance
+                0,             // slave overlap tolerance
+                true,          // Rescale weighting factors.  Bug fix, MB.
+//                 ggiInterpolation::AABB
+                overlapGgiInterpolation::BB_OCTREE  // Octree search, MB.
+
             );
 
         // Abort immediatly if uncovered faces are present
@@ -223,6 +228,15 @@ Foam::overlapGgiPolyPatch::patchToPatch() const
 
 void Foam::overlapGgiPolyPatch::calcReconFaceCellCentres() const
 {
+    if (reconFaceCellCentresPtr_)
+    {
+        FatalErrorIn
+        (
+            "void overlapGgiPolyPatch::calcReconFaceCellCentres() const"
+        )   << "Reconstructed cell centres already calculated"
+            << abort(FatalError);
+    }
+
     // Create neighbouring face centres using interpolation
     if (master())
     {
@@ -280,9 +294,15 @@ void Foam::overlapGgiPolyPatch::checkDefinition() const
 
 void Foam::overlapGgiPolyPatch::clearGeom()
 {
+    // Note: Currently deleting patch-to-patch interpolation together with
+    // expanded master and slave patches on mesh motion to avoid problems
+    // with motion of points in primitive patch.
+    // HJ, 4/Jul/2011
     deleteDemandDrivenData(expandedMasterPtr_);
     deleteDemandDrivenData(expandedSlavePtr_);
+
     deleteDemandDrivenData(patchToPatchPtr_);
+
     deleteDemandDrivenData(reconFaceCellCentresPtr_);
 }
 
@@ -294,6 +314,7 @@ void Foam::overlapGgiPolyPatch::clearOut()
     deleteDemandDrivenData(localParallelPtr_);
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::vectorField&
@@ -310,6 +331,13 @@ Foam::overlapGgiPolyPatch::reconFaceCellCentres() const
 
 void Foam::overlapGgiPolyPatch::initAddressing()
 {
+    if (active())
+    {
+        // Calculate transforms for correct GGI interpolator cut
+        calcTransforms();
+        localParallel();
+    }
+
     polyPatch::initAddressing();
 }
 
@@ -322,6 +350,17 @@ void Foam::overlapGgiPolyPatch::calcAddressing()
 
 void Foam::overlapGgiPolyPatch::initGeometry()
 {
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+    if (active())
+    {
+        // Note: Only master calculates recon; slave uses master interpolation
+        if (master())
+        {
+            reconFaceCellCentres();
+        }
+    }
+
     polyPatch::initGeometry();
 }
 
@@ -329,23 +368,29 @@ void Foam::overlapGgiPolyPatch::initGeometry()
 void Foam::overlapGgiPolyPatch::calcGeometry()
 {
     polyPatch::calcGeometry();
-
-    // Note: Calculation of transforms must be forced before the
-    // reconFaceCellCentres in order to correctly set the transformation
-    // in the interpolation routines.
-    // HJ, 3/Jul/2009
-    calcTransforms();
-
-    // Reconstruct the cell face centres
-    if (patchToPatchPtr_ && master())
-    {
-        reconFaceCellCentres();
-    }
 }
 
 
 void Foam::overlapGgiPolyPatch::initMovePoints(const pointField& p)
 {
+    clearGeom();
+
+    // Calculate transforms on mesh motion?
+    calcTransforms();
+
+    // Update interpolation for new relative position of GGI interfaces
+    // Note: currently, patches and interpolation are cleared in clearGeom()
+    // HJ. 4/Jul/2011
+//     if (patchToPatchPtr_)
+//     {
+//         patchToPatchPtr_->movePoints();
+//     }
+
+    if (active() && master())
+    {
+        reconFaceCellCentres();
+    }
+
     polyPatch::initMovePoints(p);
 }
 
@@ -353,9 +398,6 @@ void Foam::overlapGgiPolyPatch::initMovePoints(const pointField& p)
 void Foam::overlapGgiPolyPatch::movePoints(const pointField& p)
 {
     polyPatch::movePoints(p);
-
-    // Force recalculation of interpolation
-    clearGeom();
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatchTemplates.C

@@ -109,21 +109,22 @@ Foam::overlapGgiPolyPatch::interpolate(const Field<Type>& pf) const
     // Expand data
     tmp<Field<Type> > expanddata = shadow().expandData(pf);
 
-    tmp<Field<Type> > tresult;
+    tmp<Field<Type> > tresult(new Field<Type>());
+    Field<Type>& result = tresult();
 
     if (master())
     {
         // Expand slave data
-        tresult= patchToPatch().slaveToMaster(expanddata);
+        result = patchToPatch().slaveToMaster(expanddata);
     }
     else
     {
         // Expand master data
-        tresult = patchToPatch().masterToSlave(expanddata);
+        result = patchToPatch().masterToSlave(expanddata);
     }
 
     // Truncate to size
-    tresult().setSize(size());
+    result.setSize(size());
 
     return tresult;
 }

src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C

@@ -231,6 +231,9 @@ void cyclicGgiFvPatchField<Type>::initInterfaceMatrixUpdate
     const Pstream::commsTypes commsType
 ) const
 {
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+
     // Get shadow face-cells and assemble shadow field
     const unallocLabelList& sfc = cyclicGgiPatch_.shadow().faceCells();
 

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C

@@ -235,6 +235,9 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
     const Pstream::commsTypes commsType
 ) const
 {
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+
     // Get shadow face-cells and assemble shadow field
     const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
 

src/finiteVolume/fields/fvPatchFields/constraint/overlapGgi/overlapGgiFvPatchField.C

@@ -158,7 +158,19 @@ void overlapGgiFvPatchField<Type>::initEvaluate
 (
     const Pstream::commsTypes commsType
 )
-{}
+{
+    if (!this->updated())
+    {
+        this->updateCoeffs();
+    }
+
+    Field<Type>::operator=
+    (
+        this->patch().weights()*this->patchInternalField()
+      + (1.0 - this->patch().weights())*this->patchNeighbourField()
+    );
+
+}
 
 
 template<class Type>
@@ -167,39 +179,27 @@ void overlapGgiFvPatchField<Type>::evaluate
     const Pstream::commsTypes
 )
 {
-    Field<Type>::operator=
+    if (!this->updated())
-    (
+    {
-        this->patch().weights()*this->patchInternalField()
+        this->updateCoeffs();
-      + (1.0 - this->patch().weights())*this->patchNeighbourField()
+    }
-    );
 }
 
 
 template<class Type>
 void overlapGgiFvPatchField<Type>::initInterfaceMatrixUpdate
-(
-    const scalarField& psiInternal,
-    scalarField&,
-    const lduMatrix&,
-    const scalarField&,
-    const direction,
-    const Pstream::commsTypes commsType
-) const
-{}
-
-
-// Return matrix product for coupled boundary
-template<class Type>
-void overlapGgiFvPatchField<Type>::updateInterfaceMatrix
 (
     const scalarField& psiInternal,
     scalarField& result,
     const lduMatrix&,
     const scalarField& coeffs,
     const direction cmpt,
-    const Pstream::commsTypes
+    const Pstream::commsTypes commsType
 ) const
 {
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+
     // Get shadow face-cells and assemble shadow field
     const unallocLabelList& sfc = overlapGgiPatch_.shadow().faceCells();
 
@@ -222,6 +222,20 @@ void overlapGgiFvPatchField<Type>::updateInterfaceMatrix
 }
 
 
+// Return matrix product for coupled boundary
+template<class Type>
+void overlapGgiFvPatchField<Type>::updateInterfaceMatrix
+(
+    const scalarField& psiInternal,
+    scalarField& result,
+    const lduMatrix&,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes
+) const
+{}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C

@@ -41,6 +41,12 @@ Description
 namespace Foam
 {
 
+const scalar meshToMesh::cellCentreDistanceTol
+(
+    debug::tolerances("meshToMeshCellCentreDistanceTol", 1e-3)
+);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void meshToMesh::calcAddressing()
@@ -129,15 +135,18 @@ void meshToMesh::calcAddressing()
         toMesh_.cellCentres(),
         fromMesh_,
         boundaryCell,
-        oc
+        oc,
+        true
     );
 
     forAll (toMesh_.boundaryMesh(), patchi)
     {
+
         const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
 
         if (cuttingPatches_.found(toPatch.name()))
         {
+
             boundaryAddressing_[patchi].setSize(toPatch.size());
 
             cellAddresses
@@ -146,7 +155,8 @@ void meshToMesh::calcAddressing()
                 toPatch.faceCentres(),
                 fromMesh_,
                 boundaryCell,
-                oc
+                oc,
+                false
             );
         }
         else if
@@ -231,13 +241,13 @@ void meshToMesh::cellAddresses
     const pointField& points,
     const fvMesh& fromMesh,
     const List<bool>& boundaryCell,
-    const octree<octreeDataCell>& oc
+    const octree<octreeDataCell>& oc,
+    bool forceFind
 ) const
 {
 
     label nCellsOutsideAddressing = 0;
 
-
     // the implemented search method is a simple neighbour array search.
     // It starts from a cell zero, searches its neighbours and finds one
     // which is nearer to the target point than the current position.
@@ -255,6 +265,11 @@ void meshToMesh::cellAddresses
 
     forAll (points, toI)
     {
+
+        scalar localTol = cellCentreDistanceTol;
+
+        bool isBoundary = false;
+
         // pick up target position
         const vector& p = points[toI];
 
@@ -300,6 +315,7 @@ void meshToMesh::cellAddresses
             //  either way use the octree search to find it.
             if (boundaryCell[curCell])
             {
+                isBoundary = true;
                 cellAddressing_[toI] = oc.find(p);
             }
             else
@@ -356,12 +372,51 @@ void meshToMesh::cellAddresses
                 }
             }
 
-            if(cellAddressing_[toI] < 0)
+            if (cellAddressing_[toI] < 0)
             {
-
-                cellAddressing_[toI] = curCell;
-
                 nCellsOutsideAddressing++;
+
+                if (isBoundary && forceFind)
+                {
+                    // If still not found, get the closest cell within the
+                    // specified tolerance
+
+                    forAll(fromMesh.boundary(), patchi)
+                    {
+                        const fvPatch& patch = fromMesh.boundary()[patchi];
+
+                        word name = patch.name();
+
+                        label patchID =
+                            toMesh_.boundaryMesh().findPatchID(name);
+
+                        label sizePatch = 0;
+                        if (patchID > -1)
+                        {
+                            sizePatch = toMesh_.boundary()[patchID].size();
+                        }
+
+                        if
+                        (
+                            sizePatch > 0
+                        )
+                        {
+                            forAll(patch, facei)
+                            {
+                                label celli = patch.faceCells()[facei];
+
+                                const vector& centre = fromMesh.C()[celli];
+                                if (mag(points[toI] - centre) < localTol)
+                                {
+                                    localTol = mag(points[toI] - centre);
+                                    cellAddressing_[toI] = celli;
+                                }
+
+                            }
+                        }
+
+                    }
+                }
             }
         }
     }

src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C

@@ -37,7 +37,11 @@ namespace Foam
 
 defineTypeNameAndDebug(meshToMesh, 0);
 
-const scalar meshToMesh::directHitTol = 1e-5;
+const scalar meshToMesh::directHitTol
+(
+    debug::tolerances("meshToMeshdirectHitTol", 1e-5)
+);
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H

@@ -100,7 +100,8 @@ class meshToMesh
             const pointField& points,
             const fvMesh& fromMesh,
             const List<bool>& boundaryCell,
-            const octree<octreeDataCell>& oc
+            const octree<octreeDataCell>& oc,
+            bool forceFind = false
         ) const;
 
         void calculateInverseDistanceWeights() const;
@@ -113,6 +114,9 @@ class meshToMesh
         //- Direct hit tolerance
         static const scalar directHitTol;
 
+        //- Cell centre distance required to establish a hit
+        static const scalar cellCentreDistanceTol;
+
 
 public:
 
@@ -212,7 +216,7 @@ public:
             {
                 return fromMesh_;
             }
- 
+
             const fvMesh& toMesh() const
             {
                 return toMesh_;