ENH: surfaceBooleanFeatures: Update and start handling baffles

applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C

@@ -215,6 +215,29 @@ bool intersectSurfaces
 }
 
 
+label calcNormalDirection
+(
+    const vector& normal,
+    const vector& otherNormal,
+    const vector& edgeDir,
+    const vector& faceCentre,
+    const vector& pointOnEdge
+)
+{
+    vector cross = (normal ^ edgeDir);
+    cross /= mag(cross);
+
+    vector fC0tofE0 = faceCentre - pointOnEdge;
+    fC0tofE0 /= mag(fC0tofE0);
+
+    label nDir = ((cross & fC0tofE0) > 0.0 ? 1 : -1);
+
+    nDir *= ((otherNormal & fC0tofE0) > 0.0 ? -1 : 1);
+
+    return nDir;
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::noParallel();
@@ -222,6 +245,18 @@ int main(int argc, char *argv[])
     argList::validArgs.append("surface file");
     argList::validArgs.append("surface file");
 
+    argList::addBoolOption
+    (
+        "surf1Baffle",
+        "Mark surface 1 as a baffle"
+    );
+
+    argList::addBoolOption
+    (
+        "surf2Baffle",
+        "Mark surface 2 as a baffle"
+    );
+
     argList::addBoolOption
     (
         "perturb",
@@ -250,10 +285,10 @@ int main(int argc, char *argv[])
     {
         FatalErrorIn(args.executable())
             << "Unsupported action " << action << endl
-            << "Supported actions:" << validActions.toc() << exit(FatalError);
+            << "Supported actions:" << validActions.toc() << abort(FatalError);
     }
 
-    fileName surf1Name(args[2]);
+    fileName surf1Name(args.args()[2]);
     Info<< "Reading surface " << surf1Name << endl;
     triSurface surf1(surf1Name);
 
@@ -269,6 +304,9 @@ int main(int argc, char *argv[])
     surf2.writeStats(Info);
     Info<< endl;
 
+    const bool surf1Baffle = args.optionFound("surf1Baffle");
+    const bool surf2Baffle = args.optionFound("surf2Baffle");
+
     edgeIntersections edge1Cuts;
     edgeIntersections edge2Cuts;
 
@@ -330,53 +368,143 @@ int main(int argc, char *argv[])
 
     label nFeatEds = inter.cutEdges().size();
 
-    vectorField normals(2*nFeatEds, vector::zero);
+    DynamicList<vector> normals(2*nFeatEds);
     vectorField edgeDirections(nFeatEds, vector::zero);
-    labelListList edgeNormals(nFeatEds, labelList(2, label(-1)));
+    DynamicList<label> normalVolumeTypes(2*nFeatEds);
+    List<DynamicList<label> > edgeNormals(nFeatEds);
+    List<DynamicList<label> > normalDirections(nFeatEds);
 
-    triSurfaceSearch querySurf1(surf1);
+    forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
-    triSurfaceSearch querySurf2(surf2);
-
-    OFstream normalFile(sFeatFileName + "_normals.obj");
-
-    scalar scale = 0.05*min
-    (
-        querySurf1.tree().bb().mag(),
-        querySurf2.tree().bb().mag()
-    );
-
-    forAll(inter.cutEdges(), i)
     {
-        const edge& fE(inter.cutEdges()[i]);
+        const label& cutEdgeI = iter();
+        const labelPair& facePair = iter.key();
 
-        point fEC = fE.centre(inter.cutPoints());
+        const edge& fE = inter.cutEdges()[cutEdgeI];
 
-        pointIndexHit nearest1 = querySurf1.tree().findNearest(fEC, sqr(GREAT));
+        const vector& norm1 = surf1.faceNormals()[facePair.first()];
-        pointIndexHit nearest2 = querySurf2.tree().findNearest(fEC, sqr(GREAT));
+        const vector& norm2 = surf2.faceNormals()[facePair.second()];
 
-        normals[2*i] = surf1.faceNormals()[nearest1.index()];
+        DynamicList<label>& eNormals = edgeNormals[cutEdgeI];
-        normals[2*i + 1] = surf2.faceNormals()[nearest2.index()];
+        DynamicList<label>& nDirections = normalDirections[cutEdgeI];
 
-        edgeNormals[i][0] = 2*i;
+        edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
-        edgeNormals[i][1] = 2*i + 1;
 
-        edgeDirections[i] = fE.vec(inter.cutPoints());
+        normals.append(norm1);
 
+        if (surf1Baffle)
         {
-            meshTools::writeOBJ(normalFile, inter.cutPoints()[fE.start()]);
+            normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
-            meshTools::writeOBJ(normalFile, inter.cutPoints()[fE.end()]);
 
-            normalFile<< "l " << (5*i) + 1 << " " << (5*i) + 2<< endl;
+            nDirections.append(1);
+        }
+        else
+        {
+            normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+            nDirections.append
+            (
+                calcNormalDirection
+                (
+                    norm1,
+                    norm2,
+                    edgeDirections[cutEdgeI],
+                    surf1[facePair.first()].centre(surf1.points()),
+                    inter.cutPoints()[fE.start()]
+                )
+            );
+        }
 
-            meshTools::writeOBJ(normalFile, fEC);
+        eNormals.append(normals.size() - 1);
-            meshTools::writeOBJ(normalFile, fEC + scale*normals[2*i]);
-            meshTools::writeOBJ(normalFile, fEC + scale*normals[2*i + 1]);
 
-            normalFile<< "l " << (5*i) + 3 << " " << (5*i) + 4 << endl;
+        normals.append(norm2);
-            normalFile<< "l " << (5*i) + 3 << " " << (5*i) + 5 << endl;
+
+        if (surf2Baffle)
+        {
+            normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+            nDirections.append(1);
+        }
+        else
+        {
+            normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+            nDirections.append
+            (
+                calcNormalDirection
+                (
+                    norm2,
+                    norm1,
+                    edgeDirections[cutEdgeI],
+                    surf2[facePair.second()].centre(surf2.points()),
+                    inter.cutPoints()[fE.start()]
+                )
+            );
+        }
+
+        eNormals.append(normals.size() - 1);
+
+
+        if (surf1Baffle)
+        {
+            normals.append(norm2);
+
+            if (surf2Baffle)
+            {
+                normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+                nDirections.append(1);
+            }
+            else
+            {
+                normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+                nDirections.append
+                (
+                    calcNormalDirection
+                    (
+                        norm2,
+                        norm1,
+                        edgeDirections[cutEdgeI],
+                        surf2[facePair.second()].centre(surf2.points()),
+                        inter.cutPoints()[fE.start()]
+                    )
+                );
+            }
+
+            eNormals.append(normals.size() - 1);
+        }
+
+        if (surf2Baffle)
+        {
+            normals.append(norm1);
+
+            if (surf1Baffle)
+            {
+                normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+                nDirections.append(1);
+            }
+            else
+            {
+                normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+                nDirections.append
+                (
+                    calcNormalDirection
+                    (
+                        norm1,
+                        norm2,
+                        edgeDirections[cutEdgeI],
+                        surf1[facePair.first()].centre(surf1.points()),
+                        inter.cutPoints()[fE.start()]
+                    )
+                );
+            }
+
+            eNormals.append(normals.size() - 1);
         }
     }
 
+
     label internalStart = -1;
 
     if (validActions[action] == booleanSurface::UNION)
@@ -422,6 +550,19 @@ int main(int argc, char *argv[])
     // Flat, open or multiple edges are assumed to be impossible
     // Region edges are not explicitly supported by surfaceIntersection
 
+    vectorField normalsTmp(normals);
+    PackedList<2> normalVolumeTypesTmp(normalVolumeTypes);
+    labelListList edgeNormalsTmp(edgeNormals.size());
+    forAll(edgeNormalsTmp, i)
+    {
+        edgeNormalsTmp[i] = edgeNormals[i];
+    }
+    labelListList normalDirectionsTmp(normalDirections.size());
+    forAll(normalDirectionsTmp, i)
+    {
+        normalDirectionsTmp[i] = normalDirections[i];
+    }
+
     extendedFeatureEdgeMesh feMesh
     (
         IOobject
@@ -442,9 +583,11 @@ int main(int argc, char *argv[])
         nFeatEds,           // flatStart,
         nFeatEds,           // openStart,
         nFeatEds,           // multipleStart,
-        normals,
+        normalsTmp,
+        normalVolumeTypesTmp,
         edgeDirections,
-        edgeNormals,
+        normalDirectionsTmp,
+        edgeNormalsTmp,
         labelListList(0),   // featurePointNormals,
         labelListList(0),   // featurePointEdges,
         labelList(0)        // regionEdges