Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -245,6 +245,10 @@ snapControls
 
         //- Use castellatedMeshControls::features (default = true)
         explicitFeatureSnap true;
+
+        //- Detect features between multiple surfaces
+        //  (only for explicitFeatureSnap, default = false)
+        multiRegionFeatureSnap false;
 }
 
 // Settings for the layer addition.

applications/utilities/miscellaneous/foamDebugSwitches/Make/options

@@ -60,7 +60,6 @@ EXE_LIBS = \
     -lsurfaceFilmModels \
     -lsurfMesh \
     -lsystemCall \
-    -lthermalPorousZone \
     -lthermophysicalFunctions \
     -ltopoChangerFvMesh \
     -ltriSurface \

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -25,7 +25,8 @@ Application
     surfaceFeatureExtract
 
 Description
-    Extracts and writes surface features to file
+    Extracts and writes surface features to file. All but the basic feature
+    extraction is WIP.
 
 \*---------------------------------------------------------------------------*/
 
@@ -472,6 +473,176 @@ void unmarkBaffles
 }
 
 
+//- Divide into multiple normal bins
+//  - return REGION if != 2 normals
+//  - return REGION if 2 normals that make feature angle
+//  - otherwise return NONE and set normals,bins
+surfaceFeatures::edgeStatus checkFlatRegionEdge
+(
+    const triSurface& surf,
+    const scalar tol,
+    const scalar includedAngle,
+    const label edgeI
+)
+{
+    const edge& e = surf.edges()[edgeI];
+    const labelList& eFaces = surf.edgeFaces()[edgeI];
+
+    // Bin according to normal
+
+    DynamicList<Foam::vector> normals(2);
+    DynamicList<labelList> bins(2);
+
+    forAll(eFaces, eFaceI)
+    {
+        const Foam::vector& n = surf.faceNormals()[eFaces[eFaceI]];
+
+        // Find the normal in normals
+        label index = -1;
+        forAll(normals, normalI)
+        {
+            if (mag(n&normals[normalI]) > (1-tol))
+            {
+                index = normalI;
+                break;
+            }
+        }
+
+        if (index != -1)
+        {
+            bins[index].append(eFaceI);
+        }
+        else if (normals.size() >= 2)
+        {
+            // Would be third normal. Mark as feature.
+            //Pout<< "** at edge:" << surf.localPoints()[e[0]]
+            //    << surf.localPoints()[e[1]]
+            //    << " have normals:" << normals
+            //    << " and " << n << endl;
+            return surfaceFeatures::REGION;
+        }
+        else
+        {
+            normals.append(n);
+            bins.append(labelList(1, eFaceI));
+        }
+    }
+
+
+    // Check resulting number of bins
+    if (bins.size() == 1)
+    {
+        // Note: should check here whether they are two sets of faces
+        // that are planar or indeed 4 faces al coming together at an edge.
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]]
+        //    << " have single normal:" << normals[0]
+        //    << endl;
+        return surfaceFeatures::NONE;
+    }
+    else
+    {
+        // Two bins. Check if normals make an angle
+
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]] << nl
+        //    << "    normals:" << normals << nl
+        //    << "    bins   :" << bins << nl
+        //    << endl;
+
+        if (includedAngle >= 0)
+        {
+            scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+            forAll(eFaces, i)
+            {
+                const Foam::vector& ni = surf.faceNormals()[eFaces[i]];
+                for (label j=i+1; j<eFaces.size(); j++)
+                {
+                    const Foam::vector& nj = surf.faceNormals()[eFaces[j]];
+                    if (mag(ni & nj) < minCos)
+                    {
+                        //Pout<< "have sharp feature between normal:" << ni
+                        //    << " and " << nj << endl;
+
+                        // Is feature. Keep as region or convert to
+                        // feature angle? For now keep as region.
+                        return surfaceFeatures::REGION;
+                    }
+                }
+            }
+        }
+
+
+        // So now we have two normals bins but need to make sure both
+        // bins have the same regions in it.
+
+         // 1. store + or - region number depending
+        //    on orientation of triangle in bins[0]
+        const labelList& bin0 = bins[0];
+        labelList regionAndNormal(bin0.size());
+        forAll(bin0, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin0[i]]];
+            int dir = t.edgeDirection(e);
+
+            if (dir > 0)
+            {
+                regionAndNormal[i] = t.region()+1;
+            }
+            else if (dir == 0)
+            {
+                FatalErrorIn("problem.")
+                    << exit(FatalError);
+            }
+            else
+            {
+                regionAndNormal[i] = -(t.region()+1);
+            }
+        }
+
+        // 2. check against bin1
+        const labelList& bin1 = bins[1];
+        labelList regionAndNormal1(bin1.size());
+        forAll(bin1, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]];
+            int dir = t.edgeDirection(e);
+
+            label myRegionAndNormal;
+            if (dir > 0)
+            {
+                myRegionAndNormal = t.region()+1;
+            }
+            else
+            {
+                myRegionAndNormal = -(t.region()+1);
+            }
+
+            regionAndNormal1[i] = myRegionAndNormal;
+
+            label index = findIndex(regionAndNormal, -myRegionAndNormal);
+            if (index == -1)
+            {
+                // Not found.
+                //Pout<< "cannot find region " << myRegionAndNormal
+                //    << " in regions " << regionAndNormal << endl;
+
+                return surfaceFeatures::REGION;
+            }
+        }
+
+        // Passed all checks, two normal bins with the same contents.
+        //Pout<< "regionAndNormal:" << regionAndNormal << endl;
+        //Pout<< "myRegionAndNormal:" << regionAndNormal1 << endl;
+
+        return surfaceFeatures::NONE;
+    }
+}
+
+
 void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
 {
     os  << "    points : " << fem.points().size() << nl
@@ -610,6 +781,8 @@ int main(int argc, char *argv[])
 
         surfaceFeatures set(surf);
 
+        scalar includedAngle = -1;
+
         if (extractionMethod == "extractFromFile")
         {
             const fileName featureEdgeFile =
@@ -630,8 +803,7 @@ int main(int argc, char *argv[])
         }
         else if (extractionMethod == "extractFromSurface")
         {
-            const scalar includedAngle =
+            includedAngle = readScalar
-                readScalar
             (
                 surfaceDict.subDict("extractFromSurfaceCoeffs").lookup
                 (
@@ -727,13 +899,27 @@ int main(int argc, char *argv[])
             if (!nonManifoldEdges)
             {
                 Info<< "Removing all non-manifold edges"
-                    << " (edges with > 2 connected faces)" << endl;
+                    << " (edges with > 2 connected faces) unless they"
+                    << " cross multiple regions" << endl;
 
                 forAll(edgeStat, edgeI)
                 {
-                    if (surf.edgeFaces()[edgeI].size() > 2)
+                    const labelList& eFaces = surf.edgeFaces()[edgeI];
+
+                    if
+                    (
+                        eFaces.size() > 2
+                     && edgeStat[edgeI] == surfaceFeatures::REGION
+                     && (eFaces.size() % 2) == 0
+                    )
                     {
-                        edgeStat[edgeI] = surfaceFeatures::NONE;
+                        edgeStat[edgeI] = checkFlatRegionEdge
+                        (
+                            surf,
+                            1e-5,   //tol,
+                            includedAngle,
+                            edgeI
+                        );
                     }
                 }
             }

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -69,7 +69,8 @@ surface2.nas
         // Keep edges outside the box:
         outsideBox          (0 0 0)(1 1 1);
 
-        // Keep nonManifold edges (edges with >2 connected faces)
+        // Keep nonManifold edges (edges with >2 connected faces where
+        // the faces form more than two different normal planes)
         nonManifoldEdges       yes;
 
         // Keep open edges (edges with 1 connected face)

src/OSspecific/POSIX/regExp.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,12 +80,18 @@ void Foam::regExp::set(const char* pattern, const bool ignoreCase) const
             cflags |= REG_ICASE;
         }
 
-        if (regcomp(preg_, pattern, cflags) != 0)
+        int err = regcomp(preg_, pattern, cflags);
+
+        if (err != 0)
         {
+            char errbuf[200];
+            regerror(err, preg_, errbuf, sizeof(errbuf));
+
             FatalErrorIn
             (
                 "regExp::set(const char*)"
             )   << "Failed to compile regular expression '" << pattern << "'"
+                << nl << errbuf
                 << exit(FatalError);
         }
     }

src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C

@@ -266,7 +266,7 @@ Foam::Istream& Foam::ISstream::read(token& t)
                 else
                 {
                     t = sPtr;
-                    t.type() = token::STRING;
+                    t.type() = token::VARIABLE;
                 }
                 return *this;
             }

src/OpenFOAM/db/IOstreams/token/token.H

@@ -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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,6 +77,7 @@ public:
 
         PUNCTUATION,
         WORD,
+        VARIABLE,
         STRING,
         VERBATIMSTRING,
         LABEL,
@@ -331,6 +332,8 @@ public:
             inline bool isWord() const;
             inline const word& wordToken() const;
 
+            inline bool isVariable() const;
+
             inline bool isString() const;
             inline const string& stringToken() const;
 

src/OpenFOAM/db/IOstreams/token/tokenI.H

@@ -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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,7 @@ inline void token::clear()
     {
         delete wordTokenPtr_;
     }
-    else if (type_ == STRING || type_ == VERBATIMSTRING)
+    else if (type_ == STRING || type_ == VARIABLE || type_ == VERBATIMSTRING)
     {
         delete stringTokenPtr_;
     }
@@ -88,6 +88,7 @@ inline token::token(const token& t)
         break;
 
         case STRING:
+        case VARIABLE:
         case VERBATIMSTRING:
             stringTokenPtr_ = new string(*t.stringTokenPtr_);
         break;
@@ -235,14 +236,19 @@ inline const word& token::wordToken() const
     }
 }
 
+inline bool token::isVariable() const
+{
+    return (type_ == VARIABLE);
+}
+
 inline bool token::isString() const
 {
-    return (type_ == STRING || type_ == VERBATIMSTRING);
+    return (type_ == STRING || type_ == VARIABLE || type_ == VERBATIMSTRING);
 }
 
 inline const string& token::stringToken() const
 {
-    if (type_ == STRING || type_ == VERBATIMSTRING)
+    if (type_ == STRING || type_ == VARIABLE || type_ == VERBATIMSTRING)
     {
         return *stringTokenPtr_;
     }
@@ -411,6 +417,7 @@ inline void token::operator=(const token& t)
         break;
 
         case STRING:
+        case VARIABLE:
         case VERBATIMSTRING:
             stringTokenPtr_ = new string(*t.stringTokenPtr_);
         break;
@@ -518,6 +525,7 @@ inline bool token::operator==(const token& t) const
             return *wordTokenPtr_ == *t.wordTokenPtr_;
 
         case STRING:
+        case VARIABLE:
         case VERBATIMSTRING:
             return *stringTokenPtr_ == *t.stringTokenPtr_;
 
@@ -552,7 +560,11 @@ inline bool token::operator==(const word& w) const
 
 inline bool token::operator==(const string& s) const
 {
-    return ((type_ == STRING || type_ == VERBATIMSTRING) && stringToken() == s);
+    return
+    (
+        (type_ == STRING || type_ == VARIABLE || type_ == VERBATIMSTRING)
+     && stringToken() == s
+    );
 }
 
 inline bool token::operator==(const label l) const

src/OpenFOAM/db/IOstreams/token/tokenIO.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,11 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const token& t)
             os << *t.stringTokenPtr_;
         break;
 
+        case token::VARIABLE:
+            // Write variable as word so without ""
+            os.writeQuoted(*t.stringTokenPtr_, false);
+        break;
+
         case token::LABEL:
             os << t.labelToken_;
         break;
@@ -157,6 +162,10 @@ ostream& Foam::operator<<(ostream& os, const InfoProxy<token>& ip)
             os  << " the string " << t.stringToken();
         break;
 
+        case token::VARIABLE:
+            os  << " the variable " << t.stringToken();
+        break;
+
         case token::VERBATIMSTRING:
             os  << " the verbatim string " << t.stringToken();
         break;
@@ -231,6 +240,10 @@ Ostream& operator<<(Ostream& os, const InfoProxy<token>& ip)
             os  << " the string " << t.stringToken();
         break;
 
+        case token::VARIABLE:
+            os  << " the variable " << t.stringToken();
+        break;
+
         case token::VERBATIMSTRING:
             os  << " the verbatim string " << t.stringToken();
         break;

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C

@@ -55,7 +55,7 @@ void Foam::primitiveEntry::append
             newElmt(tokenIndex()++) = currToken;
         }
     }
-    else if (currToken.isString())
+    else if (currToken.isVariable())
     {
         const string& w = currToken.stringToken();
 

src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H

@@ -125,7 +125,7 @@ public:
         //- Construct from dictionary - note table is not populated
         TableBase(const word& name, const dictionary& dict);
 
-        //- Copy constructor
+        //- Copy constructor. Note: steals interpolator, tableSamples
         TableBase(const TableBase<Type>& tbl);
 
 

tutorials/mesh/snappyHexMesh/flange/system/snappyHexMeshDict

@@ -186,6 +186,10 @@ snapControls
 
         //- Use castellatedMeshControls::features
         explicitFeatureSnap true;
+
+        //- Detect features between multiple surfaces
+        //  (only for explicitFeatureSnap, default = false)
+        multiRegionFeatureSnap true;
 }