...
 
Commits (2)
  • mattijs's avatar
    ENH: faMesh: initial support for 2D topology. Unchecked. · ec114c21
    mattijs authored
    primitiveMesh:
        - mesh.nTopologicalD() : can be 3 or 2 (if all faces are two vertices)
        - if 2D:
            faceCentres : middle of edge
            faceAreas   : edge vector (currently wrongly oriented as edge;
                          should be consistent with owner/neighbour)
            cellCentres : middle of cell (using triange decomposition). Unchecked.
            cellVolumes : area of triangles. Unchecked.
    checkMesh:
        - prints 'Mesh has 2 dimensional topology'
    ec114c21
  • mattijs's avatar
    ENH: polyMesh2D: utility to collapse meshes · 17cbbd0f
    mattijs authored
    17cbbd0f
collapse2DMesh.C
EXE = $(FOAM_APPBIN)/collapse2DMesh
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
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //
......@@ -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);
......
......@@ -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_;
}
// ************************************************************************* //
......@@ -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;
......
......@@ -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);
}
}
......
......@@ -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;
}
}
}
}
......
......@@ -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,
......
......@@ -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
......