applications/utilities/mesh/generation/collapse2DMesh/Make/files

+collapse2DMesh.C
+
+EXE = $(FOAM_APPBIN)/collapse2DMesh

applications/utilities/mesh/generation/collapse2DMesh/Make/options

+EXE_INC = \
+    -DFULLDEBUG -g -O0 \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude \
+    -Iextrude2DMesh/lnInclude \
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
+
+EXE_LIBS = \
+    -lsurfMesh \
+    -ldynamicMesh \
+    -lreconstruct \
+    -lextrude2DMesh \
+    -lextrudeModel

applications/utilities/mesh/generation/collapse2DMesh/collapse2DMesh.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    collapse2DMesh
+
+Group
+    grpMeshGenerationUtilities
+
+Description
+    Takes 3D mesh with empty front and back and generate 2D version
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "polyTopoChange.H"
+#include "topoSet.H"
+#include "processorMeshes.H"
+#include "emptyPolyPatch.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Converts 2D polyMesh into polyMesh2D"
+    );
+
+    //argList::addArgument("patch", "The (empty) patch giving the geometry");
+    argList::validArgs.append("patch");
+
+    #include "addOverwriteOption.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createPolyMesh.H"
+
+    const word oldInstance = mesh.pointsInstance();
+
+    const word patchName = args[1];
+    //const bool overwrite = args.found("overwrite");
+    const bool overwrite = args.optionFound("overwrite");
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    const label keepPatchi = pbm.findPatchID(patchName);
+    if (keepPatchi == -1 || !isA<emptyPolyPatch>(pbm[keepPatchi]))
+    {
+        FatalErrorInFunction
+            << "Patch " << patchName << " does not exist or is not empty"
+            << exit(FatalError);
+    }
+
+    // Constructing a 2D mesh
+    polyTopoChange meshMod(mesh, true, 2);
+
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
+    {
+        face f(2);
+        f[0] = mesh.faces()[facei][0];
+        f[1] = mesh.faces()[facei][1];
+        meshMod.modifyFace
+        (
+            f,
+            facei,
+            own[facei],
+            nei[facei],
+            false,
+            -1,
+            -1,
+            false
+        );
+    }
+    forAll(pbm, patchi)
+    {
+        const polyPatch& pp = pbm[patchi];
+        if (!isA<emptyPolyPatch>(pp))   // || patchi == keepPatchi)
+        {
+            forAll(pp, patchFacei)
+            {
+                face f(2);
+                f[0] = pp[patchFacei][0];
+                f[1] = pp[patchFacei][1];
+                label meshFacei = pp.start()+patchFacei;
+                meshMod.modifyFace
+                (
+                    f,
+                    meshFacei,
+                    own[meshFacei],
+                    -1,
+                    false,
+                    patchi,
+                    -1,
+                    false
+                );
+            }
+        }
+        else
+        {
+            forAll(pp, patchFacei)
+            {
+                meshMod.removeFace(pp.start()+patchFacei, -1);
+            }
+        }
+    }
+
+    // Create a mesh from topo changes.
+    autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
+
+    mesh.updateMesh(morphMap);
+
+    if (!overwrite)
+    {
+        runTime++;
+    }
+    else
+    {
+        mesh.setInstance(oldInstance);
+    }
+
+    Info<< "\nWriting collapsed mesh to time = " << runTime.timeName()
+        << nl << endl;
+
+    mesh.write();
+    //topoSet::removeFiles(mesh);   // only (boundary)faceSets are wrong
+    processorMeshes::removeFiles(mesh);
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -51,6 +51,19 @@ Foam::label Foam::checkTopology
 
     Info<< "Checking topology..." << endl;
 
+    if (mesh.nTopologicalD() == 2)
+    {
+        Info<< "    Mesh has " << mesh.nTopologicalD()
+            << " dimensional topology - all faces have 2 vertices;"
+            << " cells are shell elements" << endl;
+    }
+    else if (mesh.nTopologicalD() == 1)
+    {
+        Info<< "    Mesh has " << mesh.nTopologicalD()
+            << " dimensional topology - all faces have 1 vertex;"
+            << " cells are line elements" << endl;
+    }
+
     // Check if the boundary definition is unique
     mesh.boundaryMesh().checkDefinition(true);
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -47,6 +47,7 @@ Foam::primitiveMesh::primitiveMesh()
     nInternalFaces_(0),
     nFaces_(0),
     nCells_(0),
+    nTopologicalD_(-1),
 
     cellShapesPtr_(nullptr),
     edgesPtr_(nullptr),
@@ -87,6 +88,7 @@ Foam::primitiveMesh::primitiveMesh
     nInternalFaces_(nInternalFaces),
     nFaces_(nFaces),
     nCells_(nCells),
+    nTopologicalD_(-1),
 
     cellShapesPtr_(nullptr),
     edgesPtr_(nullptr),
@@ -212,6 +214,8 @@ void Foam::primitiveMesh::reset
 {
     clearOut();
 
+    nTopologicalD_ = -1;
+
     nPoints_ = nPoints;
     nEdges_ = -1;
     nInternal0Edges_ = -1;
@@ -343,4 +347,44 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
 }
 
 
+Foam::label Foam::primitiveMesh::nTopologicalD() const
+{
+    if (nTopologicalD_ == -1)
+    {
+        const faceList& f = faces();
+
+        if (f.size())
+        {
+            nTopologicalD_ = 0;
+        }
+
+        forAll(f, facei)
+        {
+            label sz = f[facei].size();
+            if (sz > 2)
+            {
+                nTopologicalD_ = 3;
+                break;
+            }
+            else if (sz == 2)
+            {
+                nTopologicalD_ = max(nTopologicalD_, 2);
+            }
+            else if (sz == 1)
+            {
+                nTopologicalD_ = max(nTopologicalD_, 1);
+            }
+        }
+        reduce(nTopologicalD_, maxOp<label>());
+
+        if (nTopologicalD_ == -1)
+        {
+            // Zero faces. Assume 3D
+            nTopologicalD_ = 3;
+        }
+    }
+    return nTopologicalD_;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -105,6 +105,9 @@ class primitiveMesh
             //- Number of cells
             label nCells_;
 
+            //- Dimensionality of cells
+            mutable label nTopologicalD_;
+
 
         // Shapes
 
@@ -528,6 +531,10 @@ public:
 
             // Derived mesh data
 
+                //- Return the number of valid topological dimensions in the
+                //  mesh
+                label nTopologicalD() const;
+
                 //- Return cell shapes
                 const cellShapeList& cellShapes() const;
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

@@ -84,74 +84,152 @@ void Foam::primitiveMesh::makeCellCentresAndVols
     const labelList& own = faceOwner();
     const labelList& nei = faceNeighbour();
 
-    // first estimate the approximate cell centre as the average of
-    // face centres
+    if (nTopologicalD() == 2)
+    {
+        // first estimate the approximate cell centre as the average of
+        // face centres
 
-    vectorField cEst(nCells(), Zero);
-    labelField nCellFaces(nCells(), 0);
+        vectorField cEst(nCells(), Zero);
+        labelField nCellFaces(nCells(), 0);
 
-    forAll(own, facei)
-    {
-        cEst[own[facei]] += fCtrs[facei];
-        nCellFaces[own[facei]] += 1;
-    }
+        forAll(own, facei)
+        {
+            cEst[own[facei]] += fCtrs[facei];
+            nCellFaces[own[facei]] += 1;
+        }
 
-    forAll(nei, facei)
-    {
-        cEst[nei[facei]] += fCtrs[facei];
-        nCellFaces[nei[facei]] += 1;
-    }
+        forAll(nei, facei)
+        {
+            cEst[nei[facei]] += fCtrs[facei];
+            nCellFaces[nei[facei]] += 1;
+        }
 
-    forAll(cEst, celli)
-    {
-        cEst[celli] /= nCellFaces[celli];
-    }
+        forAll(cEst, celli)
+        {
+            cEst[celli] /= nCellFaces[celli];
+        }
 
-    forAll(own, facei)
-    {
-        // Calculate 3*face-pyramid volume
-        scalar pyr3Vol =
-            fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
+        forAll(own, facei)
+        {
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol = mag
+            (
+                fAreas[facei] ^ (fCtrs[facei] - cEst[own[facei]])
+            );
 
-        // Calculate face-pyramid centre
-        vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
+            // Calculate face-pyramid centre
+            vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
 
-        // Accumulate volume-weighted face-pyramid centre
-        cellCtrs[own[facei]] += pyr3Vol*pc;
+            // Accumulate volume-weighted face-pyramid centre
+            cellCtrs[own[facei]] += pyr3Vol*pc;
 
-        // Accumulate face-pyramid volume
-        cellVols[own[facei]] += pyr3Vol;
-    }
+            // Accumulate face-pyramid volume
+            cellVols[own[facei]] += pyr3Vol;
+        }
 
-    forAll(nei, facei)
-    {
-        // Calculate 3*face-pyramid volume
-        scalar pyr3Vol =
-            fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
+        forAll(nei, facei)
+        {
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol = mag
+            (
+                fAreas[facei] ^ (cEst[nei[facei]] - fCtrs[facei])
+            );
 
-        // Calculate face-pyramid centre
-        vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
+            // Calculate face-pyramid centre
+            vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
 
-        // Accumulate volume-weighted face-pyramid centre
-        cellCtrs[nei[facei]] += pyr3Vol*pc;
+            // Accumulate volume-weighted face-pyramid centre
+            cellCtrs[nei[facei]] += pyr3Vol*pc;
 
-        // Accumulate face-pyramid volume
-        cellVols[nei[facei]] += pyr3Vol;
-    }
+            // Accumulate face-pyramid volume
+            cellVols[nei[facei]] += pyr3Vol;
+        }
+
+        forAll(cellCtrs, celli)
+        {
+            if (mag(cellVols[celli]) > VSMALL)
+            {
+                cellCtrs[celli] /= cellVols[celli];
+            }
+            else
+            {
+                cellCtrs[celli] = cEst[celli];
+            }
+        }
 
-    forAll(cellCtrs, celli)
+        cellVols *= (1.0/3.0);
+    }
+    else
     {
-        if (mag(cellVols[celli]) > VSMALL)
+        // first estimate the approximate cell centre as the average of
+        // face centres
+
+        vectorField cEst(nCells(), Zero);
+        labelField nCellFaces(nCells(), 0);
+
+        forAll(own, facei)
         {
-            cellCtrs[celli] /= cellVols[celli];
+            cEst[own[facei]] += fCtrs[facei];
+            nCellFaces[own[facei]] += 1;
         }
-        else
+
+        forAll(nei, facei)
         {
-            cellCtrs[celli] = cEst[celli];
+            cEst[nei[facei]] += fCtrs[facei];
+            nCellFaces[nei[facei]] += 1;
         }
-    }
 
-    cellVols *= (1.0/3.0);
+        forAll(cEst, celli)
+        {
+            cEst[celli] /= nCellFaces[celli];
+        }
+
+        forAll(own, facei)
+        {
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol =
+                fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
+
+            // Calculate face-pyramid centre
+            vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
+
+            // Accumulate volume-weighted face-pyramid centre
+            cellCtrs[own[facei]] += pyr3Vol*pc;
+
+            // Accumulate face-pyramid volume
+            cellVols[own[facei]] += pyr3Vol;
+        }
+
+        forAll(nei, facei)
+        {
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol =
+                fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
+
+            // Calculate face-pyramid centre
+            vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
+
+            // Accumulate volume-weighted face-pyramid centre
+            cellCtrs[nei[facei]] += pyr3Vol*pc;
+
+            // Accumulate face-pyramid volume
+            cellVols[nei[facei]] += pyr3Vol;
+        }
+
+        forAll(cellCtrs, celli)
+        {
+            if (mag(cellVols[celli]) > VSMALL)
+            {
+                cellCtrs[celli] /= cellVols[celli];
+            }
+            else
+            {
+                cellCtrs[celli] = cEst[celli];
+            }
+        }
+
+        cellVols *= (1.0/3.0);
+    }
 }
 
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C

@@ -79,56 +79,75 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
 {
     const faceList& fs = faces();
 
-    forAll(fs, facei)
+    if (nTopologicalD() == 2)
     {
-        const labelList& f = fs[facei];
-        label nPoints = f.size();
-
-        // If the face is a triangle, do a direct calculation for efficiency
-        // and to avoid round-off error-related problems
-        if (nPoints == 3)
-        {
-            fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
-            fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
-        }
-        else
+        forAll(fs, facei)
         {
-            vector sumN = Zero;
-            scalar sumA = 0.0;
-            vector sumAc = Zero;
+            const labelList& f = fs[facei];
 
-            point fCentre = p[f[0]];
-            for (label pi = 1; pi < nPoints; pi++)
+            if (f.size() != 2)
             {
-                fCentre += p[f[pi]];
-            }
-
-            fCentre /= nPoints;
-
-            for (label pi = 0; pi < nPoints; pi++)
-            {
-                const point& nextPoint = p[f[(pi + 1) % nPoints]];
-
-                vector c = p[f[pi]] + nextPoint + fCentre;
-                vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
-                scalar a = mag(n);
-
-                sumN += n;
-                sumA += a;
-                sumAc += a*c;
+                FatalErrorInFunction
+                    << "Illegal face " << f << " at "
+                    << UIndirectList<point>(p, f) << endl;
             }
+            fCtrs[facei] = 0.5*(p[f[0]] + p[f[1]]);
+            fAreas[facei] = p[f[1]]-p[f[0]];
+        }
+    }
+    else
+    {
+        forAll(fs, facei)
+        {
+            const labelList& f = fs[facei];
+            label nPoints = f.size();
 
-            // This is to deal with zero-area faces. Mark very small faces
-            // to be detected in e.g., processorPolyPatch.
-            if (sumA < ROOTVSMALL)
+            // If the face is a triangle, do a direct calculation for efficiency
+            // and to avoid round-off error-related problems
+            if (nPoints == 3)
             {
-                fCtrs[facei] = fCentre;
-                fAreas[facei] = Zero;
+                fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
+                fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
             }
             else
             {
-                fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
-                fAreas[facei] = 0.5*sumN;
+                vector sumN = Zero;
+                scalar sumA = 0.0;
+                vector sumAc = Zero;
+
+                point fCentre = p[f[0]];
+                for (label pi = 1; pi < nPoints; pi++)
+                {
+                    fCentre += p[f[pi]];
+                }
+
+                fCentre /= nPoints;
+
+                for (label pi = 0; pi < nPoints; pi++)
+                {
+                    const point& nextPoint = p[f[(pi + 1) % nPoints]];
+
+                    vector c = p[f[pi]] + nextPoint + fCentre;
+                    vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
+                    scalar a = mag(n);
+
+                    sumN += n;
+                    sumA += a;
+                    sumAc += a*c;
+                }
+
+                // This is to deal with zero-area faces. Mark very small faces
+                // to be detected in e.g., processorPolyPatch.
+                if (sumA < ROOTVSMALL)
+                {
+                    fCtrs[facei] = fCentre;
+                    fAreas[facei] = Zero;
+                }
+                else
+                {
+                    fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
+                    fAreas[facei] = 0.5*sumN;
+                }
             }
         }
     }

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -404,7 +404,14 @@ void Foam::polyTopoChange::checkFace
         }
     }
 
-    if (f.size() < 3 || findIndex(f, -1) != -1)
+    if
+    (
+        (
+            (nTopologicalD_ == 3 && f.size() < 3)
+         || (nTopologicalD_ == 2 && f.size() != 2)
+        )
+     || findIndex(f, -1) != -1
+    )
     {
         FatalErrorInFunction
             << "Illegal vertices in face"
@@ -1100,7 +1107,14 @@ void Foam::polyTopoChange::compact
             //labelList oldF(f);
             renumberCompact(localPointMap, f);
 
-            if (!faceRemoved(facei) && f.size() < 3)
+            if
+            (
+               !faceRemoved(facei)
+             && (
+                    (nTopologicalD_ == 3 && f.size() < 3)
+                 || (nTopologicalD_ == 2 && f.size() != 2)
+                )
+            )
             {
                 FatalErrorInFunction
                     << "Created illegal face " << f
@@ -2181,9 +2195,15 @@ void Foam::polyTopoChange::compactAndReorder
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
+Foam::polyTopoChange::polyTopoChange
+(
+    const label nPatches,
+    const bool strict,
+    const label nTopologicalD
+)
 :
     strict_(strict),
+    nTopologicalD_(nTopologicalD),
     nPatches_(nPatches),
     points_(0),
     pointMap_(0),
@@ -2208,17 +2228,24 @@ Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
     cellFromEdge_(0),
     cellFromFace_(0),
     cellZone_(0)
-{}
+{
+    if (nTopologicalD_ != 2 && nTopologicalD_ != 3)
+    {
+        FatalErrorInFunction << "Only 3D and 2D supported" << exit(FatalError);
+    }
+}
 
 
 // Construct from components
 Foam::polyTopoChange::polyTopoChange
 (
     const polyMesh& mesh,
-    const bool strict
+    const bool strict,
+    const label nTopologicalD
 )
 :
     strict_(strict),
+    nTopologicalD_(nTopologicalD),
     nPatches_(0),
     points_(0),
     pointMap_(0),
@@ -2244,6 +2271,10 @@ Foam::polyTopoChange::polyTopoChange
     cellFromFace_(0),
     cellZone_(0)
 {
+    if (nTopologicalD_ != 2 && nTopologicalD_ != 3)
+    {
+        FatalErrorInFunction << "Only 3D and 2D supported" << exit(FatalError);
+    }
     addMesh
     (
         mesh,

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

@@ -103,6 +103,9 @@ class polyTopoChange
         //  when adding/removing data.
         bool strict_;
 
+        //- Type of mesh to generate: 3D (default) or 2D (faMesh)
+        label nTopologicalD_;
+
 
         // Patches
 
@@ -408,10 +411,20 @@ public:
 
         //- Construct without mesh. Either specify nPatches or use
         //  setNumPatches before trying to make a mesh (makeMesh, changeMesh)
-        polyTopoChange(const label nPatches, const bool strict = true);
+        polyTopoChange
+        (
+            const label nPatches,
+            const bool strict = true,
+            const label nTopologicalD = 3
+        );
 
         //- Construct from mesh. Adds all points/face/cells from mesh.
-        polyTopoChange(const polyMesh& mesh, const bool strict = true);
+        polyTopoChange
+        (
+            const polyMesh& mesh,
+            const bool strict = true,
+            const label nTopologicalD = 3
+        );
 
 
     // Member Functions