Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/conversion/ansysToFoam/ansysToFoam.L

@@ -63,11 +63,16 @@ label maxNodei = 0;
 
 SLPtrList<labelList> slCellLabels;
 SLList<label> slCellMap;
+SLList<label> slCellType;
 label maxCelli = 0;
 
 PtrList<SLList<label> > slPatchCells;
 PtrList<SLList<label> > slPatchCellFaces;
 
+// Cell types
+Map<word> cellTypes;
+label currentTypei = -1;
+
 
 // Dummy yywrap to keep yylex happy at compile time.
 // It is called by yylex but is not used as the mechanism to change file.
@@ -108,6 +113,8 @@ value                 {floatNum}
 node                  ^{space}"N"{cspace}
 element               ^{space}"EN"{cspace}
 bface                 ^{space}"SFE"{cspace}
+elementTypeName       ^{space}"ET"{cspace}
+elementType           ^{space}"TYPE"{cspace}
 
 
 %%
@@ -160,6 +167,7 @@ bface                 ^{space}"SFE"{cspace}
 
         slCellMap.append(celli);
         slCellLabels.append(new labelList(labels));
+        slCellType.append(currentTypei);
     }
 
 
@@ -210,6 +218,37 @@ bface                 ^{space}"SFE"{cspace}
     }
 
 
+{elementTypeName}{label}{cspace}{identifier}{space}\n {
+
+        IStringStream elementStream(YYText());
+        char tag,c;
+        label cellTypei;
+        word cellTypeName;
+        elementStream
+            >> tag >> tag           // skip 'ET'
+            >> c >> cellTypei
+            >> c >> cellTypeName;
+
+        Info<< "Read typeName " << cellTypeName
+            << " for type " << cellTypei << endl;
+
+        cellTypes.insert(cellTypei, cellTypeName);
+    }
+
+
+{elementType}{label}{space}\n {
+        IStringStream elementStream(YYText());
+        char tag,c;
+        label cellTypei;
+        elementStream
+            >> tag >> tag >> tag >> tag     // skip 'TYPE'
+            >> c >> cellTypei;
+
+        currentTypei = cellTypei;
+    }
+
+
+
  /* ------------------------------------------------------------------------- *\
     ------ Ignore remaining space and \n s.  Any other characters are errors.
  \* ------------------------------------------------------------------------- */
@@ -231,6 +270,29 @@ bface                 ^{space}"SFE"{cspace}
 #include <fstream>
 using std::ifstream;
 
+
+label findFace(const polyMesh& mesh, const face& f)
+{
+    const labelList& pFaces = mesh.pointFaces()[f[0]];
+
+    forAll(pFaces, i)
+    {
+        label faceI = pFaces[i];
+
+        if (mesh.faces()[faceI] == f)
+        {
+            return faceI;
+        }
+    }
+
+    FatalErrorIn("findFace(const polyMesh&, const face&)")
+        << "Cannot find a face matching " << f
+        << exit(FatalError);
+
+    return -1;
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::noParallel();
@@ -253,7 +315,7 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 
-    const fileName ansysFile = args[1];
+    fileName ansysFile(args.additionalArgs()[0]);
     ifstream ansysStream(ansysFile.c_str());
 
     if (!ansysStream)
@@ -377,6 +439,34 @@ int main(int argc, char *argv[])
         }
     }
 
+
+    const word defaultFacesName = "defaultFaces";
+    word defaultFacesType = emptyPolyPatch::typeName;
+
+    // Create dummy mesh just to find out what are internal/external
+    // faces
+    autoPtr<polyMesh> dummyMesh
+    (
+        new polyMesh
+        (
+            IOobject
+            (
+                "dummyMesh",
+                runTime.constant(),
+                runTime
+            ),
+            xferCopy(points),
+            cellShapes,
+            faceListList(0),
+            wordList(0),
+            wordList(0),
+            defaultFacesName,
+            defaultFacesType,
+            wordList(0)
+        )
+    );
+
+
     // Warning: tet face order has changed between version 1.9.6 and 2.0
     //
     label faceIndex[7][6] =
@@ -423,10 +513,53 @@ int main(int argc, char *argv[])
 
         boundary[patchI] = patchFaces;
         patchNames[patchI] = word("patch") + name(patchI + 1);
+    }
+
+
+    //
+    // Lookup the face labels for all the boundary faces
+    //
+    labelListList boundaryFaceLabels(boundary.size());
+    forAll(boundary, patchI)
+    {
+        const faceList& bFaces = boundary[patchI];
+        labelList& bFaceLabels = boundaryFaceLabels[patchI];
+        bFaceLabels.setSize(bFaces.size());
+        forAll(bFaces, i)
+        {
+            bFaceLabels[i] = findFace(dummyMesh(), bFaces[i]);
+        }
+    }
+
+
+    // Now split the boundary faces into external and internal faces. All
+    // faces go into faceZones and external faces go into patches.
+    List<faceList> patchFaces(slPatchCells.size());
+    labelList patchNFaces(slPatchCells.size(), 0);
+    forAll(boundary, patchI)
+    {
+        const faceList& bFaces = boundary[patchI];
+        const labelList& bFaceLabels = boundaryFaceLabels[patchI];
+
+        patchFaces[patchI].setSize(bFaces.size());
+
+        forAll(bFaces, i)
+        {
+            if (!dummyMesh().isInternalFace(bFaceLabels[i]))
+            {
+                patchFaces[patchI][patchNFaces[patchI]++] = bFaces[i];
+            }
+        }
+        patchFaces[patchI].setSize(patchNFaces[patchI]);
+
         Info<< "Patch " << patchI << " named " << patchNames[patchI]
             << ": " << boundary[patchI].size() << " faces" << endl;
     }
 
+    // We no longer need the dummyMesh
+    dummyMesh.clear();
+
+
     Info<< "ansysToFoam: " << endl
         << "Ansys file format does not provide information about the type of "
         << "the patch (eg. wall, symmetry plane, cyclic etc)." << endl
@@ -434,10 +567,8 @@ int main(int argc, char *argv[])
         << "as type patch. Please reset after mesh conversion as necessary."
         << endl;
 
-    wordList patchTypes(boundary.size(), polyPatch::typeName);
-    word defaultFacesName = "defaultFaces";
-    word defaultFacesType = emptyPolyPatch::typeName;
-    wordList patchPhysicalTypes(boundary.size());
+    wordList patchTypes(patchFaces.size(), polyPatch::typeName);
+    wordList patchPhysicalTypes(patchFaces.size());
 
     preservePatchTypes
     (
@@ -461,7 +592,7 @@ int main(int argc, char *argv[])
         ),
         xferMove(points),
         cellShapes,
-        boundary,
+        patchFaces,
         patchNames,
         patchTypes,
         defaultFacesName,
@@ -469,6 +600,90 @@ int main(int argc, char *argv[])
         patchPhysicalTypes
     );
 
+
+    if (cellTypes.size() > 0 || patchNames.size() > 0)
+    {
+        DynamicList<pointZone*> pz;
+        DynamicList<faceZone*> fz;
+        DynamicList<cellZone*> cz;
+
+        // FaceZones
+        forAll(boundaryFaceLabels, patchI)
+        {
+            if (boundaryFaceLabels[patchI].size())
+            {
+                // Re-do the boundaryFaceLabels since the boundary face
+                // labels will be different on the pShapeMesh.
+                const faceList& bFaces = boundary[patchI];
+                labelList& bFaceLabels = boundaryFaceLabels[patchI];
+                forAll(bFaceLabels, i)
+                {
+                    bFaceLabels[i] = findFace(pShapeMesh, bFaces[i]);
+                }
+
+                Info<< "Creating faceZone " <<  patchNames[patchI]
+                    << " with " << bFaceLabels.size() << " faces" << endl;
+
+                fz.append
+                (
+                    new faceZone
+                    (
+                        patchNames[patchI],
+                        bFaceLabels,
+                        boolList(bFaceLabels.size(), false),
+                        fz.size(),
+                        pShapeMesh.faceZones()
+                    )
+                );
+            }
+        }
+
+
+        // CellZones
+        labelList types = cellTypes.sortedToc();
+
+        forAll(types, j)
+        {
+            label cellType = types[j];
+
+            // Pick up cells in zone
+            DynamicList<label> addr;
+
+            SLList<label>::iterator cellMapIter = slCellMap.begin();
+            SLList<label>::iterator typeIter = slCellType.begin();
+
+            for
+            (
+                ;
+                typeIter != slCellType.end();
+                ++typeIter, ++cellMapIter
+            )
+            {
+                if (typeIter() == cellType)
+                {
+                    addr.append(cellMap[cellMapIter()]);
+                }
+            }
+
+            Info<< "Creating cellZone " << cellTypes[cellType]
+                << " with " << addr.size() << " cells" << endl;
+
+            cz.append
+            (
+                new cellZone
+                (
+                    cellTypes[cellType],
+                    addr,
+                    j,
+                    pShapeMesh.cellZones()
+                )
+            );
+        }
+
+        pShapeMesh.addZones(pz, fz, cz);
+    }
+
+
     // Set the precision of the points data to 10
     IOstream::defaultPrecision(10);
 

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -7,32 +7,6 @@
 #include "pointSet.H"
 #include "IOmanip.H"
 
-bool Foam::checkSync(const wordList& names)
-{
-    List<wordList> allNames(Pstream::nProcs());
-    allNames[Pstream::myProcNo()] = names;
-    Pstream::gatherList(allNames);
-    Pstream::scatterList(allNames);
-
-    bool hasError = false;
-
-    for (label procI = 1; procI < allNames.size(); procI++)
-    {
-        if (allNames[procI] != allNames[0])
-        {
-            hasError = true;
-
-            Info<< " ***Inconsistent zones across processors, "
-                   "processor 0 has zones:" << allNames[0]
-                << ", processor " << procI << " has zones:"
-                << allNames[procI]
-                << endl;
-        }
-    }
-    return hasError;
-}
-
-
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -51,35 +25,22 @@ Foam::label Foam::checkTopology
     mesh.boundaryMesh().checkParallelSync(true);
 
     // Check names of zones are equal
-    if (checkSync(mesh.cellZones().names()))
+    mesh.cellZones().checkDefinition(true);
+    if (mesh.cellZones().checkParallelSync(true))
     {
         noFailedChecks++;
     }
-    if (checkSync(mesh.faceZones().names()))
+    mesh.faceZones().checkDefinition(true);
+    if (mesh.faceZones().checkParallelSync(true))
     {
         noFailedChecks++;
     }
-    if (checkSync(mesh.pointZones().names()))
+    mesh.pointZones().checkDefinition(true);
+    if (mesh.pointZones().checkParallelSync(true))
     {
         noFailedChecks++;
     }
 
-    // Check contents of faceZones consistent
-    {
-        forAll(mesh.faceZones(), zoneI)
-        {
-            if (mesh.faceZones()[zoneI].checkParallelSync(false))
-            {
-                Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
-                    << " is not correctly synchronised"
-                    << " across coupled boundaries."
-                    << " (coupled faces are either not both "
-                    << " present in set or have same flipmap)" << endl;
-                noFailedChecks++;
-            }
-        }
-    }
-
     {
         pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
         if (mesh.checkPoints(true, &points))

applications/utilities/mesh/manipulation/checkMesh/checkTopology.H

@@ -5,7 +5,5 @@ namespace Foam
 {
     class polyMesh;
 
-    bool checkSync(const wordList& names);
-
     label checkTopology(const polyMesh&, const bool, const bool);
 }

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -60,6 +60,8 @@ void Foam::ODESolver::solve
 
     scalar x = xStart;
     scalar h = hEst;
+    scalar hNext = 0;
+    scalar hPrev = 0;
 
     for (label nStep=0; nStep<MAXSTP; nStep++)
     {
@@ -73,14 +75,24 @@ void Foam::ODESolver::solve
         if ((x + h - xEnd)*(x + h - xStart) > 0.0)
         {
             h = xEnd - x;
+            hPrev = hNext;
         }
 
-        scalar hNext, hDid;
+        hNext = 0;
+        scalar hDid;
         solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext);
 
         if ((x - xEnd)*(xEnd - xStart) >= 0.0)
         {
-            hEst = hNext;
+            if (hPrev != 0)
+            {
+                hEst = hPrev;
+            }
+            else
+            {
+                hEst = hNext;
+            }
+
             return;
         }
 

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C

@@ -27,6 +27,7 @@ License
 #include "entry.H"
 #include "demandDrivenData.H"
 #include "stringListOps.H"
+#include "Pstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -372,6 +373,84 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkDefinition
 }
 
 
+template<class ZoneType, class MeshType>
+bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
+(
+    const bool report
+) const
+{
+    if (!Pstream::parRun())
+    {
+        return false;
+    }
+
+
+    const PtrList<ZoneType>& zones = *this;
+
+    bool hasError = false;
+
+    // Collect all names
+    List<wordList> allNames(Pstream::nProcs());
+    allNames[Pstream::myProcNo()] = this->names();
+    Pstream::gatherList(allNames);
+    Pstream::scatterList(allNames);
+
+    List<wordList> allTypes(Pstream::nProcs());
+    allTypes[Pstream::myProcNo()] = this->types();
+    Pstream::gatherList(allTypes);
+    Pstream::scatterList(allTypes);
+
+    // Have every processor check but only master print error.
+
+    for (label procI = 1; procI < allNames.size(); procI++)
+    {
+        if
+        (
+            (allNames[procI] != allNames[0])
+         || (allTypes[procI] != allTypes[0])
+        )
+        {
+            hasError = true;
+
+            if (debug || (report && Pstream::master()))
+            {
+                Info<< " ***Inconsistent zones across processors, "
+                       "processor 0 has zone names:" << allNames[0]
+                    << " zone types:" << allTypes[0]
+                    << " processor " << procI << " has zone names:"
+                    << allNames[procI]
+                    << " zone types:" << allTypes[procI]
+                    << endl;
+            }
+        }
+    }
+
+    // Check contents
+    if (!hasError)
+    {
+        forAll(zones, zoneI)
+        {
+            if (zones[zoneI].checkParallelSync(false))
+            {
+                hasError = true;
+
+                if (debug || (report && Pstream::master()))
+                {
+                    Info<< " ***Zone " << zones[zoneI].name()
+                        << " of type " << zones[zoneI].type()
+                        << " is not correctly synchronised"
+                        << " across coupled boundaries."
+                        << " (coupled faces are either not both "
+                        << " present in set or have same flipmap)" << endl;
+                }
+            }
+        }
+    }
+
+    return hasError;
+}
+
+
 // Correct zone mesh after moving points
 template<class ZoneType, class MeshType>
 void Foam::ZoneMesh<ZoneType, MeshType>::movePoints(const pointField& p)

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H

@@ -149,6 +149,10 @@ public:
         //- Check zone definition. Return true if in error.
         bool checkDefinition(const bool report = false) const;
 
+        //- Check whether all procs have all zones and in same order. Return
+        //  true if in error.
+        bool checkParallelSync(const bool report = false) const;
+
         //- Correct zone mesh after moving points
         void movePoints(const pointField&);
 

src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H

@@ -209,6 +209,13 @@ public:
         //- Check zone definition. Return true if in error.
         virtual bool checkDefinition(const bool report = false) const;
 
+        //- Check whether zone is synchronised across coupled boundaries. Return
+        //  true if in error.
+        virtual bool checkParallelSync(const bool report = false) const
+        {
+            return false;
+        }
+
         //- Write dictionary
         virtual void writeDict(Ostream&) const;
 

src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.C

@@ -29,6 +29,7 @@ License
 #include "polyMesh.H"
 #include "primitiveMesh.H"
 #include "demandDrivenData.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -134,6 +135,49 @@ bool Foam::pointZone::checkDefinition(const bool report) const
 }
 
 
+bool Foam::pointZone::checkParallelSync(const bool report) const
+{
+    const polyMesh& mesh = zoneMesh().mesh();
+
+    labelList maxZone(mesh.nPoints(), -1);
+    labelList minZone(mesh.nPoints(), labelMax);
+    forAll(*this, i)
+    {
+        label pointI = operator[](i);
+        maxZone[pointI] = index();
+        minZone[pointI] = index();
+    }
+    syncTools::syncPointList(mesh, maxZone, maxEqOp<label>(), -1);
+    syncTools::syncPointList(mesh, minZone, minEqOp<label>(), labelMax);
+
+    bool error = false;
+
+    forAll(maxZone, pointI)
+    {
+        // Check point in zone on both sides
+        if (maxZone[pointI] != minZone[pointI])
+        {
+            if (report && !error)
+            {
+                Info<< " ***Problem with pointZone " << index()
+                    << " named " << name()
+                    << ". Point " << pointI
+                    << " at " << mesh.points()[pointI]
+                    << " is in zone "
+                    << (minZone[pointI] == labelMax ? -1 : minZone[pointI])
+                    << " on some processors and in zone "
+                    << maxZone[pointI]
+                    << " on some other processors."
+                    << endl;
+            }
+            error = true;
+        }
+    }
+
+    return error;
+}
+
+
 void Foam::pointZone::writeDict(Ostream& os) const
 {
     os  << nl << name_ << nl << token::BEGIN_BLOCK << nl

src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H

@@ -208,6 +208,10 @@ public:
         //- Check zone definition. Return true if in error.
         virtual bool checkDefinition(const bool report = false) const;
 
+        //- Check whether zone is synchronised across coupled boundaries. Return
+        //  true if in error.
+        virtual bool checkParallelSync(const bool report = false) const;
+