Merge remote-tracking branch 'origin/nr/HenrikRusche' into nextRelease

src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/addCellLayer.C

@@ -52,7 +52,7 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     if (debug)
     {
-        Pout<< "void layerAdditionRemoval::addCellLayer("
+        Info<< "void layerAdditionRemoval::addCellLayer("
             << "polyTopoChange& ref) const for object " << name() << " : "
             << "Adding cell layer" << endl;
     }
@@ -74,7 +74,7 @@ void Foam::layerAdditionRemoval::addCellLayer
     {
         if (debug)
         {
-            Pout<< "void layerAdditionRemoval::addCellLayer("
+            Info<< "void layerAdditionRemoval::addCellLayer("
                 << "polyTopoChange& ref) const "
                 << " for object " << name() << " : "
                 << "Using edges for point insertion" << endl;
@@ -93,7 +93,7 @@ void Foam::layerAdditionRemoval::addCellLayer
     {
         if (debug)
         {
-            Pout<< "void layerAdditionRemoval::addCellLayer("
+            Info<< "void layerAdditionRemoval::addCellLayer("
                 << "polyTopoChange& ref) const "
                 << " for object " << name() << " : "
                 << "A valid layer could not be found in front of "
@@ -550,7 +550,7 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     if (debug)
     {
-        Pout<< "void layerAdditionRemoval::addCellLayer("
+        Info<< "void layerAdditionRemoval::addCellLayer("
             << "polyTopoChange& ref) const "
             << " for object " << name() << " : "
             << "Finished adding cell layer" << endl;

src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C

@@ -328,7 +328,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
 
     if (debug)
     {
-        Pout<< "bool layerAdditionRemoval::changeTopology() const "
+        Info<< "bool layerAdditionRemoval::changeTopology() const "
             << " for object " << name() << " : " << nl
             << "Layer thickness: min: " << minDelta
             << " max: " << maxDelta << " avg: " << avgDelta
@@ -345,7 +345,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     {
         if (debug)
         {
-            Pout << "First step. No addition/removal" << endl;
+            Info<< "First step. No addition/removal" << endl;
         }
 
         // No topological changes allowed before first mesh motion
@@ -372,7 +372,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
 
                     if (debug)
                     {
-                        Pout<< "bool layerAdditionRemoval::changeTopology() "
+                        Info<< "bool layerAdditionRemoval::changeTopology() "
                             << " const for object " << name() << " : "
                             << "Triggering layer removal" << endl;
                     }
@@ -404,7 +404,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
         {
             if (debug)
             {
-                Pout<< "bool layerAdditionRemoval::changeTopology() const "
+                Info<< "bool layerAdditionRemoval::changeTopology() const "
                     << " for object " << name() << " : "
                     << "Triggering layer addition" << endl;
             }
@@ -439,7 +439,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Info<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer removal. " << endl;
         }
@@ -455,7 +455,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Info<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer addition. " << endl;
         }
@@ -470,17 +470,17 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
 {
     if (debug)
     {
-        Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
+        Info<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
             << " for object " << name() << " : "
             << "Clearing addressing on external request. ";
 
         if (pointsPairingPtr_ || facesPairingPtr_)
         {
-            Pout << "Pointers set." << endl;
+            Info << "Pointers set." << endl;
         }
         else
         {
-            Pout << "Pointers not set." << endl;
+            Info << "Pointers not set." << endl;
         }
     }
 

src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/setLayerPairing.C

@@ -123,11 +123,15 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
 //     << "curLocalFace: " << curLocalFace << nl
 //     << "lidFace: " << lidFace
 //     << " master index: " << lidFace.masterIndex()
-//     << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
+//     << " oppositeIndex: " << lidFace.oppositeIndex() << nl
+//     << "faces[oppositeIndex]: " << faces[lidFace.oppositeIndex()] << endl;
 
         // Grab the opposite face for face collapse addressing
         ftc[faceI] = lidFace.oppositeIndex();
 
+        // lidFace is only valid in simple cases
+        const face& lidFace2 = faces[lidFace.oppositeIndex()];
+
         // Using the local face insert the points into the lid list
         forAll (curLocalFace, pointI)
         {
@@ -149,25 +153,25 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
                     // arbitray polyhedra!
                     // HR, 1/May/2011
 
-                    label curPoint = faces[mf[faceI]][pointI];
+                    const label curPoint = faces[mf[faceI]][pointI];
 
                     const labelListList& ppAddr = mesh.pointPoints();
 
-                    const labelList& p1List = ppAddr[lidFace[pointI]];
+                    const labelList& p1List = ppAddr[lidFace2[pointI]];
                     const labelList& p2List = ppAddr[ptc[clp]];
 
                     bool found = false;
-                    forAll(p1List, p1I)
+                    forAll (p1List, p1I)
                     {
-                        label p1 = p1List[p1I];
+                        const label p1 = p1List[p1I];
 
                         if (p1 != curPoint)
                         {
-                            forAll(p2List, p2I)
+                            forAll (p2List, p2I)
                             {
                                 label p2 = p2List[p2I];
 
-                                if(p1 == p2)
+                                if (p1 == p2)
                                 {
                                     ptc[clp] = p1;
                                     found = true;
@@ -176,7 +180,10 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
                             }
                         }
 
-                        if(found) { break; }
+                        if (found)
+                        {
+                            break;
+                        }
                     }
 
                     if (!found)

src/foam/interpolations/GGIInterpolation/GGIInterpolation.C

@@ -202,8 +202,17 @@ GGIInterpolation<MasterPatch, SlavePatch>::uncoveredSlaveFaces() const
 
 
 template<class MasterPatch, class SlavePatch>
-bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints()
+bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints
+(
+    const tensorField& forwardT,
+    const tensorField& reverseT,
+    const vectorField& forwardSep
+)
 {
+    this->forwardT_ = forwardT;
+    this->reverseT_ = reverseT;
+    this->forwardSep_ = forwardSep;
+
     clearOut();
 
     return true;

src/foam/interpolations/GGIInterpolation/GGIInterpolation.H

@@ -133,18 +133,18 @@ class GGIInterpolation
         //  master plane.  Size equals number of slave faces; zero length
         //  indicates no transform.  Size 1 indicates constant transform
         //  HJ, 6/Jan/2009
-        const tensorField forwardT_;
+        tensorField forwardT_;
 
         //- Master-to-slave transformation tensor.  Transforms slave data to
         //  master plane.  Size equals number of master faces; zero length
         //  indicates no transform.  Size 1 indicates constant transform
         //  HJ, 6/Jan/2009
-        const tensorField reverseT_;
+        tensorField reverseT_;
 
         //- Slave-to-master separation vector. Translation of slave data to
         //  master plane.  Size equals number of slave faces; zero length
         //  indicates no translation.  MB, 28/Jan/2009
-        const vectorField forwardSep_;
+        vectorField forwardSep_;
 
         //- Master non-overlap face tolerance factor
         const scalar masterNonOverlapFaceTol_;
@@ -552,7 +552,12 @@ public:
     // Edit
 
         //- Correct weighting factors for moving mesh.
-        bool movePoints();
+        bool movePoints
+        (
+            const tensorField& forwardT,
+            const tensorField& reverseT,
+            const vectorField& forwardSep
+        );
 };
 
 

src/foam/meshes/meshShapes/cell/oppositeCellFace.C

@@ -68,7 +68,7 @@ Foam::label Foam::cell::opposingFaceLabel
         {
             bool sharedPoint = false;
 
-            // Compare every vertex of the current face agains the
+            // Compare every vertex of the current face against the
             // vertices of the master face
             forAll (curFace, pointI)
             {
@@ -103,6 +103,7 @@ Foam::label Foam::cell::opposingFaceLabel
                         Info<< "Multiple faces not sharing vertex: "
                             << oppositeFaceLabel << " and "
                             << curFaceLabels[faceI] << endl;
+
                         return -1;
                     }
                 }
@@ -179,6 +180,9 @@ Foam::oppositeFace Foam::cell::opposingFace
                                 break;
                             }
                         }
+
+                        // None of them matches, return a failure
+                        return oppositeFace(face(0), masterFaceLabel, oppFaceLabel);
                     }
                 }
             }

src/foam/meshes/meshShapes/face/oppositeFace.H

@@ -95,7 +95,7 @@ public:
         //- Does the opposite face exist?
         inline bool found() const
         {
-            return oppositeIndex_ >= 0;
+            return size() > 0;
         }
 };
 

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

@@ -833,7 +833,12 @@ void Foam::ggiPolyPatch::initMovePoints(const pointField& p)
     // Update interpolation for new relative position of GGI interfaces
     if (patchToPatchPtr_)
     {
-        patchToPatchPtr_->movePoints();
+        patchToPatchPtr_->movePoints
+        (
+            forwardT(),
+            reverseT(),
+            -separation()
+        );
     }
 
     // Recalculate send and receive maps

src/foam/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -199,11 +199,13 @@ void Foam::processorPolyPatch::calcGeometry()
             scalar nbrMagSf = mag(neighbFaceAreas_[facei]);
             scalar avSf = (magSf + nbrMagSf)/2.0;
 
-            if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
+            if (magSf < SMALL && nbrMagSf < SMALL)
             {
                 // Undetermined normal. Use dummy normal to force separation
-                // check. (note use of sqrt(VSMALL) since that is how mag
+                // check.
-                // scales)
+                // (note use of sqrt(VSMALL) since that is how mag scales)
+                // HR, 11/12/2013: Found a face with area = 1e-21 before
+                // topo deflation. Hence must use SMALL here.
                 faceNormals[facei] = point(1, 0, 0);
                 nbrFaceNormals[facei] = faceNormals[facei];
             }

src/foam/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C

@@ -1001,7 +1001,12 @@ void Foam::regionCouplePolyPatch::initMovePoints(const pointField& p)
     // Update interpolation for new relative position of GGI interfaces
     if (patchToPatchPtr_)
     {
-        patchToPatchPtr_->movePoints();
+        patchToPatchPtr_->movePoints
+        (
+            forwardT(),
+            reverseT(),
+            shadow().separation()
+        );
     }
 
     // Recalculate send and receive maps

src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C

@@ -413,7 +413,7 @@ solidContactFvPatchVectorField::solidContactFvPatchVectorField
                         *masterFaceZonePatchPtr_, // to zone
                         alg_,
                         dir_
-                        );
+                    );
             }
             else if (interpolationMethod_ == "ggi")
             {
@@ -439,7 +439,7 @@ solidContactFvPatchVectorField::solidContactFvPatchVectorField
                         0.0,
                         true,
                         ggiInterpolation::AABB
-                        );
+                    );
             }
             else
             {
@@ -706,7 +706,12 @@ void solidContactFvPatchVectorField::updateCoeffs()
                 {
                     if (slaveToMasterGgiInterpolatorPtr_)
                     {
-                        slaveToMasterGgiInterpolatorPtr_->movePoints();
+                        slaveToMasterGgiInterpolatorPtr_->movePoints
+                        (
+                            tensorField(0),
+                            tensorField(0),
+                            vectorField(0)
+                        );
                     }
                     else if (slaveToMasterPatchToPatchInterpolatorPtr_)
                     {
@@ -897,14 +902,14 @@ tmp<vectorField> solidContactFvPatchVectorField::interpolateSlaveToMaster
       (
           masterFaceZonePatchPtr_->size(),
           vector::zero
-          );
+      );
   if (slaveToMasterPatchToPatchInterpolatorPtr_)
   {
       globalMasterInterpField =
           slaveToMasterPatchToPatchInterpolatorPtr_->faceInterpolate<vector>
           (
               globalSlaveField
-              );
+          );
   }
   else if (slaveToMasterGgiInterpolatorPtr_)
   {
@@ -912,7 +917,7 @@ tmp<vectorField> solidContactFvPatchVectorField::interpolateSlaveToMaster
           slaveToMasterGgiInterpolatorPtr_->slaveToMaster
           (
               globalSlaveField
-              );
+          );
   }
   else
   {
@@ -1112,7 +1117,7 @@ void solidContactFvPatchVectorField::moveFaceZonePatches()
               *masterFaceZonePatchPtr_, // to zone
               alg_,
               dir_
-              );
+          );
   }
   else if (slaveToMasterGgiInterpolatorPtr_)
   {
@@ -1132,7 +1137,7 @@ void solidContactFvPatchVectorField::moveFaceZonePatches()
               0.0,
               true,
               ggiInterpolation::AABB
-              );
+          );
   }
 
   // and primitive patch interpolators