From ec197599c4b6be2c620575c143062e1c457a2554 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 8 Aug 2008 12:12:49 +0100
Subject: [PATCH] handle internal marked faces

---
 .../mesh/conversion/tetgenToFoam/tetgenToFoam.C    | 203 +++++++++++++--------
 1 file changed, 129 insertions(+), 74 deletions(-)

diff --git a/applications/utilities/mesh/conversion/tetgenToFoam/tetgenToFoam.C b/applications/utilities/mesh/conversion/tetgenToFoam/tetgenToFoam.C
index 79f0fd3..1e50cce 100644
--- a/applications/utilities/mesh/conversion/tetgenToFoam/tetgenToFoam.C
+++ b/applications/utilities/mesh/conversion/tetgenToFoam/tetgenToFoam.C
@@ -54,8 +54,9 @@ Description
 NOTE:
 - for some reason boundary faces point inwards. I just reverse them
 always. Might use some geometric check instead.
-- marked faces might not actually be boundary faces of mesh. This is not handled
-and you'll have to run without face file (-noFaceFile option)
+- marked faces might not actually be boundary faces of mesh.
+This is hopefully handled now by first constructing without boundaries
+and then reconstructing with boundary faces.
 
 \*---------------------------------------------------------------------------*/
 
@@ -69,20 +70,40 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Find label of face.
+label findFace(const primitiveMesh& 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 primitiveMesh&, const face&)")
+        << "Cannot find face " << f << " in mesh." << abort(FatalError);
+
+    return -1;
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
     argList::validArgs.append("file prefix");
     argList::validOptions.insert("noFaceFile", "");
-    argList::validOptions.insert("overwrite", "");
 
 #   include "setRootCase.H"
 #   include "createTime.H"
 
 
     bool readFaceFile = !args.options().found("noFaceFile");
-    bool overwrite = args.options().found("overwrite");
 
     fileName prefix(args.additionalArgs()[0]);
 
@@ -287,16 +308,44 @@ int main(int argc, char *argv[])
     }
 
 
-    label nPatches = 0;
+    //
+    // Construct mesh with default boundary only
+    //
+
+    autoPtr<polyMesh> meshPtr
+    (
+        new polyMesh
+        (
+            IOobject
+            (
+                polyMesh::defaultRegion,
+                runTime.constant(),
+                runTime
+            ),
+            points,
+            cells,
+            faceListList(0),
+            wordList(0),    //boundaryPatchNames
+            wordList(0),    //boundaryPatchTypes
+            "defaultFaces",
+            polyPatch::typeName,
+            wordList(0)
+        )
+    );
+    const polyMesh& mesh = meshPtr;
 
-    // List of Foam vertices per boundary face
-    faceList boundaryFaces;
 
-    // For each boundary faces the Foam patchID
-    labelList boundaryPatch;
 
     if (readFaceFile)
     {
+        label nPatches = 0;
+
+        // List of Foam vertices per boundary face
+        faceList boundaryFaces;
+
+        // For each boundary faces the Foam patchID
+        labelList boundaryPatch;
+
         //
         // read boundary faces
         //
@@ -366,48 +415,59 @@ int main(int argc, char *argv[])
                     f[2-i] = nodeToPoint[nodeI];
                 }
 
-                boundaryFaces[faceI] = f;
 
-                if (nFaceAttr > 0)
+                if (findFace(mesh, f) >= mesh.nInternalFaces())
                 {
-                    // First attribute is the region number
-                    faceLine >> region;
+                    boundaryFaces[faceI] = f;
 
+                    if (nFaceAttr > 0)
+                    {
+                        // First attribute is the region number
+                        faceLine >> region;
 
-                    // Get Foam patchID and update region->patch table.
-                    label patchI = 0;
 
-                    Map<label>::iterator patchFind = regionToPatch.find(region);
+                        // Get Foam patchID and update region->patch table.
+                        label patchI = 0;
 
-                    if (patchFind == regionToPatch.end())
-                    {
-                        patchI = nPatches;
+                        Map<label>::iterator patchFind =
+                            regionToPatch.find(region);
 
-                        Info<< "Mapping tetgen region " << region
-                            << " to Foam patch "
-                            << patchI << endl;
+                        if (patchFind == regionToPatch.end())
+                        {
+                            patchI = nPatches;
 
-                        regionToPatch.insert(region, nPatches++);
-                    }
-                    else
-                    {
-                        patchI = patchFind();
-                    }
+                            Info<< "Mapping tetgen region " << region
+                                << " to Foam patch "
+                                << patchI << endl;
 
-                    boundaryPatch[faceI] = patchI;
+                            regionToPatch.insert(region, nPatches++);
+                        }
+                        else
+                        {
+                            patchI = patchFind();
+                        }
 
-                    // Skip remaining attributes
-                    for (label i = 1; i < nFaceAttr; i++)
-                    {
-                        faceLine >> dummy;
+                        boundaryPatch[faceI] = patchI;
+
+                        // Skip remaining attributes
+                        for (label i = 1; i < nFaceAttr; i++)
+                        {
+                            faceLine >> dummy;
+                        }
                     }
-                }
 
-                faceI++;
+                    faceI++;
+                }
             }
         }
 
-        // Print region to patch mapping
+   
+        // Trim
+        boundaryFaces.setSize(faceI);
+        boundaryPatch.setSize(faceI);
+
+
+         // Print region to patch mapping
         Info<< "Regions:" << endl;
 
         for
@@ -421,28 +481,23 @@ int main(int argc, char *argv[])
                 << iter() << endl;
         }
         Info<< endl;
-    }
 
 
-    // Storage for boundary faces
+        // Storage for boundary faces
+        faceListList patchFaces(nPatches);
+        wordList patchNames(nPatches);
 
-    faceListList patchFaces(nPatches);
-
-    wordList patchNames(nPatches);
-
-    forAll(patchNames, patchI)
-    {
-        patchNames[patchI] = word("patch") + name(patchI);
-    }
+        forAll(patchNames, patchI)
+        {
+            patchNames[patchI] = word("patch") + name(patchI);
+        }
 
-    wordList patchTypes(nPatches, polyPatch::typeName);
-    word defaultFacesName = "defaultFaces";
-    word defaultFacesType = polyPatch::typeName;
-    wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
+        wordList patchTypes(nPatches, polyPatch::typeName);
+        word defaultFacesName = "defaultFaces";
+        word defaultFacesType = polyPatch::typeName;
+        wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
 
 
-    if (readFaceFile)
-    {
         // Sort boundaryFaces by patch using boundaryPatch.
         List<DynamicList<face> > allPatchFaces(nPatches);
 
@@ -464,34 +519,34 @@ int main(int argc, char *argv[])
         }
 
         Info<< endl;
-    }
 
-    if (!overwrite)
-    {
-        runTime++;
-    }
 
-    polyMesh mesh
-    (
-        IOobject
+        meshPtr.reset
         (
-            polyMesh::defaultRegion,
-            runTime.constant(),
-            runTime
-        ),
-        points,
-        cells,
-        patchFaces,
-        patchNames,
-        patchTypes,
-        defaultFacesName,
-        defaultFacesType,
-        patchPhysicalTypes
-    );
+            new polyMesh
+            (
+                IOobject
+                (
+                    polyMesh::defaultRegion,
+                    runTime.constant(),
+                    runTime
+                ),
+                points,
+                cells,
+                patchFaces,
+                patchNames,
+                patchTypes,
+                defaultFacesName,
+                defaultFacesType,
+                patchPhysicalTypes
+            )
+        );
+    }
+
 
     Info<< "Writing mesh to " << runTime.constant() << endl << endl;
 
-    mesh.write();
+    meshPtr().write();
 
 
     Info<< "End\n" << endl;
-- 
2.11.4.GIT