From 01737e14e7ed6e4f55411ba15a7a022f143a0fc5 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 17 Oct 2018 09:13:35 +0100
Subject: [PATCH] BUG: removeFaces: do not remove if inbetween same cells. See
 #998.

This is the additional fix to make dynamic unrefinement work again.
---
 .../polyTopoChange/hexRef8/hexRef8.C          | 20 ++++++++++++-
 .../polyTopoChange/removeFaces.C              | 30 +++++++++++++++----
 2 files changed, 44 insertions(+), 6 deletions(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index 0c0a5e2480..a66bd7c486 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -5572,9 +5572,27 @@ void Foam::hexRef8::setUnrefinement
 
         if (facesToRemove.size() != splitFaces.size())
         {
+            // Dump current mesh
+            {
+                const_cast<polyMesh&>(mesh_).setInstance
+                (
+                    mesh_.time().timeName()
+                );
+                mesh_.write();
+                pointSet pSet(mesh_, "splitPoints", splitPointLabels);
+                pSet.write();
+                faceSet fSet(mesh_, "splitFaces", splitFaces);
+                fSet.write();
+                faceSet removeSet(mesh_, "facesToRemove", facesToRemove);
+                removeSet.write();
+            }
+
             FatalErrorInFunction
                 << "Ininitial set of split points to unrefine does not"
                 << " seem to be consistent or not mid points of refined cells"
+                << " splitPoints:" << splitPointLabels.size()
+                << " splitFaces:" << splitFaces.size()
+                << " facesToRemove:" << facesToRemove.size()
                 << abort(FatalError);
         }
     }
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
index 7baa85a3d1..496188b2ad 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
@@ -910,7 +910,7 @@ void Foam::removeFaces::setRefinement
         {
             if (nFacesPerEdge[edgeI] == 2)
             {
-                // See if they are two boundary faces
+                // Get the two face labels
                 label f0 = -1;
                 label f1 = -1;
 
@@ -920,7 +920,7 @@ void Foam::removeFaces::setRefinement
                 {
                     label facei = eFaces[i];
 
-                    if (!removedFace[facei] && !mesh_.isInternalFace(facei))
+                    if (!removedFace[facei])
                     {
                         if (f0 == -1)
                         {
@@ -934,7 +934,7 @@ void Foam::removeFaces::setRefinement
                     }
                 }
 
-                if (f0 != -1 && f1 != -1)
+                if (!mesh_.isInternalFace(f0) && !mesh_.isInternalFace(f1))
                 {
                     // Edge has two boundary faces remaining.
                     // See if should be merged.
@@ -979,7 +979,7 @@ void Foam::removeFaces::setRefinement
                         }
                     }
                 }
-                else if (f0 != -1 || f1 != -1)
+                else if (mesh_.isInternalFace(f0) != mesh_.isInternalFace(f1))
                 {
                     const edge& e = mesh_.edges()[edgeI];
 
@@ -1000,7 +1000,27 @@ void Foam::removeFaces::setRefinement
                 else
                 {
                     // Both kept faces are internal. Mark edge for preserving
-                    nFacesPerEdge[edgeI] = 3;
+                    // if inbetween different cells. If inbetween same cell
+                    // pair we probably want to merge them to
+                    //  - avoid upper-triangular ordering problems
+                    //  - allow hex unrefinement (expects single face inbetween
+                    //    cells)
+
+                    const edge ownEdge
+                    (
+                        cellRegion[mesh_.faceOwner()[f0]],
+                        cellRegion[mesh_.faceNeighbour()[f0]]
+                    );
+                    const edge neiEdge
+                    (
+                        cellRegion[mesh_.faceOwner()[f1]],
+                        cellRegion[mesh_.faceNeighbour()[f1]]
+                    );
+
+                    if (ownEdge != neiEdge)
+                    {
+                        nFacesPerEdge[edgeI] = 3;
+                    }
                 }
             }
         }