From 1f7fb08060f2fc3de05391e8b84c2ff73516701b Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Thu, 23 Nov 2023 14:55:58 +0200
Subject: [PATCH] BUG: wrong numbering when determining the zoneIDs

of the STL written by topology optimisation.

BUG: when determining which mesh faces are cut by iso-surface faces,
only append the latter if it contains more than two points
---
 .../topOVariablesBase/topOVariablesBase.C     | 49 ++++++++++---------
 .../topOVariablesBase/topOVariablesBase.H     | 10 ++--
 2 files changed, 30 insertions(+), 29 deletions(-)

diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
index d620cb7061..0299e96951 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
@@ -141,8 +141,7 @@ bool Foam::topOVariablesBase::addCuttingFaceToIsoline
     {
         // Check whether any of points of the new iso-surface face are already
         // present in the surface. To reduce the number of comparisons, only
-        // points on iso-surface faces belonging to neighbouring cells are
-        // checked
+        // points on iso-surface faces cutting the current cell are checked
         labelList uniquePointIDs(facePoints.size(), -1);
         DynamicList<point> uniqueFacePoints(facePoints.size());
         DynamicList<label> uniqueFacePointEdges(facePoints.size());
@@ -362,8 +361,8 @@ void Foam::topOVariablesBase::writeSurfaceFiles
         }
 
         // Faces passed in previous zones
-        labelList cumulZoneSizes(surfZoneIds.size(), 0);
-        for (label zi = 1; zi < surfZoneIds.size() - 1; ++zi)
+        labelList cumulZoneSizes(nSerialPatches + 1, 0);
+        for (label zi = 1; zi < cumulZoneSizes.size(); ++zi)
         {
             cumulZoneSizes[zi] = cumulZoneSizes[zi - 1] + zoneSizes[zi - 1];
         }
@@ -426,8 +425,6 @@ void Foam::topOVariablesBase::writeSurfaceFiles
 }
 
 
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::topOVariablesBase::topOVariablesBase
@@ -453,7 +450,7 @@ Foam::topOVariablesBase::topOVariablesBase
     isoSurfFolder_
         (mesh.time().globalPath()/"optimisation"/"topOIsoSurfaces"),
     meshFaceToChangedFace_(),
-    changedFacesPerCuttingFace_(),
+    //changedFacesPerCuttingFace_(),
     surfPoints_(),
     surfFaces_()
 {
@@ -615,23 +612,27 @@ void Foam::topOVariablesBase::writeFluidSolidInterface
                 }
             }
 
-            addCuttingFaceToIsoline
+            if
             (
-                cutCell.facePoints(),
-                nSerialPatches,
-                cellCutFaces,
-                cuttingFacesPerMeshFace,
-                isoSurfPts,
-                isoSurfFaces,
-                zoneIDs
-            );
-
-            for (const label facei : cellCutFaces)
-            {
-                cuttingFacesPerMeshFace[facei].push_back
+                addCuttingFaceToIsoline
                 (
-                    isoSurfFaces.size() - 1
-                );
+                    cutCell.facePoints(),
+                    nSerialPatches,
+                    cellCutFaces,
+                    cuttingFacesPerMeshFace,
+                    isoSurfPts,
+                    isoSurfFaces,
+                    zoneIDs
+                )
+            )
+            {
+                for (const label facei : cellCutFaces)
+                {
+                    cuttingFacesPerMeshFace[facei].push_back
+                    (
+                        isoSurfFaces.size() - 1
+                    );
+                }
             }
         }
     }
@@ -664,8 +665,8 @@ void Foam::topOVariablesBase::writeFluidSolidInterface
     writeSurfaceFiles(surfPoints_, surfFaces_, zoneIds, nSerialPatches);
 
     // Invert changedFace-to-cuttingFace map for the sensitivity computations
-    changedFacesPerCuttingFace_ =
-        invertOneToMany(surfFaces_.size(), changedFaceToCutFace);
+    //changedFacesPerCuttingFace_ =
+    //    invertOneToMany(surfFaces_.size(), changedFaceToCutFace);
 
     // Transform origin cut faces to a global numbering
     labelList cuttingFacesPerProc(Pstream::nProcs(), Zero);
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.H
index ee47c1b0f5..a5ac4e1779 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.H
@@ -77,7 +77,7 @@ protected:
         Map<label> meshFaceToChangedFace_;
 
         //- Per cutting face, the changed faces owned by it
-        labelListList changedFacesPerCuttingFace_;
+        //labelListList changedFacesPerCuttingFace_;
 
         //- Iso-surface points
         pointField surfPoints_;
@@ -208,10 +208,10 @@ public:
             return zones_.getBetaMax();
         }
 
-        inline const labelListList& changedFacesPerCuttingFace() const
-        {
-            return changedFacesPerCuttingFace_;
-        }
+       //inline const labelListList& changedFacesPerCuttingFace() const
+       //{
+       //    return changedFacesPerCuttingFace_;
+       //}
 
         inline const pointField& surfacePoints() const
         {