ENH: respect face orientation when decomposing polyhedra.

ENH: use face::trianglesQuads() method for PV3FoamReader as well.
- this avoids missing faces (and weird cells) in the decomposed polyhedra.
This commit is contained in:
Mark Olesen
2010-03-04 11:20:20 +01:00
parent 5726d59f63
commit 89615f708e
31 changed files with 188 additions and 157 deletions

View File

@ -159,13 +159,6 @@ Note
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
static const label VTK_TETRA = 10;
static const label VTK_PYRAMID = 14;
static const label VTK_WEDGE = 13;
static const label VTK_HEXAHEDRON = 12;
template<class GeoField>
void print(const char* msg, Ostream& os, const PtrList<GeoField>& flds)
{
@ -229,10 +222,7 @@ labelList getSelectedPatches
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
@ -285,11 +275,11 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
bool doWriteInternal = !args.optionFound("noInternal");
bool doFaceZones = !args.optionFound("noFaceZones");
bool doLinks = !args.optionFound("noLinks");
bool binary = !args.optionFound("ascii");
bool useTimeName = args.optionFound("useTimeName");
const bool doWriteInternal = !args.optionFound("noInternal");
const bool doFaceZones = !args.optionFound("noFaceZones");
const bool doLinks = !args.optionFound("noLinks");
const bool binary = !args.optionFound("ascii");
const bool useTimeName = args.optionFound("useTimeName");
if (binary && (sizeof(floatScalar) != 4 || sizeof(label) != 4))
{
@ -299,7 +289,7 @@ int main(int argc, char *argv[])
<< exit(FatalError);
}
bool nearCellValue = args.optionFound("nearCellValue");
const bool nearCellValue = args.optionFound("nearCellValue");
if (nearCellValue)
{
@ -308,7 +298,7 @@ int main(int argc, char *argv[])
<< nl << endl;
}
bool noPointValues = args.optionFound("noPointValues");
const bool noPointValues = args.optionFound("noPointValues");
if (noPointValues)
{
@ -316,7 +306,7 @@ int main(int argc, char *argv[])
<< "Outputting cell values only" << nl << endl;
}
bool allPatches = args.optionFound("allPatches");
const bool allPatches = args.optionFound("allPatches");
List<wordRe> excludePatches;
if (args.optionFound("excludePatches"))
@ -392,15 +382,8 @@ int main(int argc, char *argv[])
Info<< "Time: " << runTime.timeName() << endl;
word timeDesc = "";
if (useTimeName)
{
timeDesc = runTime.timeName();
}
else
{
timeDesc = name(runTime.timeIndex());
}
word timeDesc =
useTimeName ? runTime.timeName() : Foam::name(runTime.timeIndex());
// Check for new polyMesh/ and update mesh, fvMeshSubset and cell
// decomposition.
@ -470,10 +453,7 @@ int main(int argc, char *argv[])
IOobjectList objects(mesh, runTime.timeName());
HashSet<word> selectedFields;
if (args.optionFound("fields"))
{
args.optionLookup("fields")() >> selectedFields;
}
args.optionReadIfPresent("fields", selectedFields);
// Construct the vol fields (on the original mesh if subsetted)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -33,10 +33,10 @@ Description
#include "polyMesh.H"
#include "cellShape.H"
#include "cellModeller.H"
#include "Swap.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
:
mesh_(mesh),
@ -61,6 +61,9 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
// Number of additional cells generated by the decomposition of polyhedra
label nAddCells = 0;
// face owner is needed to determine the face orientation
const labelList& owner = mesh.faceOwner();
// Scan for cells which need to be decomposed and count additional points
// and cells
@ -71,7 +74,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
if
(
model != hex
// && model != wedge // See above.
// && model != wedge // See above.
&& model != prism
&& model != pyr
&& model != tet
@ -111,7 +114,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
cellTypes_.setSize(cellShapes.size() + nAddCells);
// Set counters for additional points and additional cells
label api = 0, aci = 0;
label addPointI = 0, addCellI = 0;
forAll(cellShapes, cellI)
{
@ -134,15 +137,13 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
}
else if (cellModel == prism)
{
vtkVerts.setSize(6);
vtkVerts[0] = cellShape[0];
vtkVerts[1] = cellShape[2];
vtkVerts[2] = cellShape[1];
vtkVerts[3] = cellShape[3];
vtkVerts[4] = cellShape[5];
vtkVerts[5] = cellShape[4];
// VTK has a different node order for VTK_WEDGE
// their triangles point outwards!
vtkVerts = cellShape;
Foam::Swap(vtkVerts[1], vtkVerts[2]);
Foam::Swap(vtkVerts[4], vtkVerts[5]);
// VTK calls this a wedge.
cellTypes_[cellI] = VTK_WEDGE;
}
else if (cellModel == tetWedge)
@ -175,34 +176,28 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
// }
else if (cellModel == hex)
{
vtkVerts.setSize(8);
vtkVerts[0] = cellShape[0];
vtkVerts[1] = cellShape[1];
vtkVerts[2] = cellShape[2];
vtkVerts[3] = cellShape[3];
vtkVerts[4] = cellShape[4];
vtkVerts[5] = cellShape[5];
vtkVerts[6] = cellShape[6];
vtkVerts[7] = cellShape[7];
vtkVerts = cellShape;
cellTypes_[cellI] = VTK_HEXAHEDRON;
}
else
{
// Polyhedral cell. Decompose into tets + prisms.
// (see dxFoamExec/createDxConnections.C)
// Mapping from additional point to cell
addPointCellLabels_[api] = cellI;
addPointCellLabels_[addPointI] = cellI;
// The new vertex from the cell-centre
const label newVertexLabel = mesh_.nPoints() + addPointI;
// Whether to insert cell in place of original or not.
bool substituteCell = true;
const labelList& cFaces = mesh_.cells()[cellI];
forAll(cFaces, cFaceI)
{
const face& f = mesh_.faces()[cFaces[cFaceI]];
const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
// Number of triangles and quads in decomposition
label nTris = 0;
@ -216,75 +211,95 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
label quadi = 0;
f.trianglesQuads(mesh_.points(), trii, quadi, triFcs, quadFcs);
forAll(quadFcs, quadi)
forAll(quadFcs, quadI)
{
label thisCellI = -1;
label thisCellI;
if (substituteCell)
{
thisCellI = cellI;
substituteCell = false;
}
else
{
thisCellI = mesh_.nCells() + aci;
superCells_[aci] = cellI;
aci++;
thisCellI = mesh_.nCells() + addCellI;
superCells_[addCellI++] = cellI;
}
labelList& addVtkVerts = vertLabels_[thisCellI];
addVtkVerts.setSize(5);
const face& quad = quadFcs[quadi];
const face& quad = quadFcs[quadI];
addVtkVerts[0] = quad[0];
addVtkVerts[1] = quad[1];
addVtkVerts[2] = quad[2];
addVtkVerts[3] = quad[3];
addVtkVerts[4] = mesh_.nPoints() + api;
// Ensure we have the correct orientation for the
// base of the primitive cell shape.
// If the cell is face owner, the orientation needs to be
// flipped.
// At the moment, VTK doesn't actually seem to care if
// negative cells are defined, but we'll do it anyhow
// (for safety).
if (isOwner)
{
addVtkVerts[0] = quad[3];
addVtkVerts[1] = quad[2];
addVtkVerts[2] = quad[1];
addVtkVerts[3] = quad[0];
}
else
{
addVtkVerts[0] = quad[0];
addVtkVerts[1] = quad[1];
addVtkVerts[2] = quad[2];
addVtkVerts[3] = quad[3];
}
addVtkVerts[4] = newVertexLabel;
cellTypes_[thisCellI] = VTK_PYRAMID;
}
forAll(triFcs, trii)
forAll(triFcs, triI)
{
label thisCellI = -1;
label thisCellI;
if (substituteCell)
{
thisCellI = cellI;
substituteCell = false;
}
else
{
thisCellI = mesh_.nCells() + aci;
superCells_[aci] = cellI;
aci++;
thisCellI = mesh_.nCells() + addCellI;
superCells_[addCellI++] = cellI;
}
labelList& addVtkVerts = vertLabels_[thisCellI];
const face& tri = triFcs[trii];
const face& tri = triFcs[triI];
addVtkVerts.setSize(4);
addVtkVerts[0] = tri[0];
addVtkVerts[1] = tri[1];
addVtkVerts[2] = tri[2];
addVtkVerts[3] = mesh_.nPoints() + api;
// See note above about the orientation.
if (isOwner)
{
addVtkVerts[0] = tri[2];
addVtkVerts[1] = tri[1];
addVtkVerts[2] = tri[0];
}
else
{
addVtkVerts[0] = tri[0];
addVtkVerts[1] = tri[1];
addVtkVerts[2] = tri[2];
}
addVtkVerts[3] = newVertexLabel;
cellTypes_[thisCellI] = VTK_TETRA;
}
}
api++;
addPointI++;
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -80,7 +80,7 @@ public:
// Public static data
// this must be consistent with the enumeration in "vtkCell.H"
//- equivalent to enumeration in "vtkCellType.h"
enum vtkTypes
{
VTK_TRIANGLE = 5,
@ -88,9 +88,10 @@ public:
VTK_QUAD = 9,
VTK_TETRA = 10,
VTK_PYRAMID = 14,
VTK_WEDGE = 13,
VTK_HEXAHEDRON = 12,
VTK_WEDGE = 13,
VTK_PYRAMID = 14,
VTK_POLYHEDRON = 42
};
// Constructors

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -91,7 +91,8 @@
animateable="0">
<BooleanDomain name="bool"/>
<Documentation>
Use vtkPolyhedron instead of decomposing polyhedra
Use vtkPolyhedron instead of decomposing polyhedra.
!!Actually uses vtkConvexPointSet until this is properly supported in VTK!!
</Documentation>
</IntVectorProperty>

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -31,6 +31,7 @@ License
#include "fvMesh.H"
#include "cellModeller.H"
#include "vtkOpenFOAMPoints.H"
#include "Swap.H"
// VTK includes
#include "vtkCellArray.h"
@ -45,6 +46,13 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
polyDecomp& decompInfo
)
{
const cellModel& tet = *(cellModeller::lookup("tet"));
const cellModel& pyr = *(cellModeller::lookup("pyr"));
const cellModel& prism = *(cellModeller::lookup("prism"));
const cellModel& wedge = *(cellModeller::lookup("wedge"));
const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
const cellModel& hex = *(cellModeller::lookup("hex"));
vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
if (debug)
@ -53,36 +61,27 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
printMemory();
}
const cellShapeList& cellShapes = mesh.cellShapes();
// Number of additional points needed by the decomposition of polyhedra
label nAddPoints = 0;
// Number of additional cells generated by the decomposition of polyhedra
label nAddCells = 0;
// face owner is needed to determine the face orientation
const labelList& owner = mesh.faceOwner();
labelList& superCells = decompInfo.superCells();
labelList& addPointCellLabels = decompInfo.addPointCellLabels();
const cellModel& tet = *(cellModeller::lookup("tet"));
const cellModel& pyr = *(cellModeller::lookup("pyr"));
const cellModel& prism = *(cellModeller::lookup("prism"));
const cellModel& wedge = *(cellModeller::lookup("wedge"));
const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
const cellModel& hex = *(cellModeller::lookup("hex"));
// Scan for cells which need to be decomposed and count additional points
// and cells
if (debug)
{
Info<< "... building cell-shapes" << endl;
}
const cellShapeList& cellShapes = mesh.cellShapes();
if (debug)
{
Info<< "... scanning" << endl;
}
// count number of cells to decompose
// Scan for cells which need to be decomposed and count additional points
// and cells
if (!reader_->GetUseVTKPolyhedron())
{
forAll(cellShapes, cellI)
@ -105,10 +104,10 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
{
const face& f = mesh.faces()[cFaces[cFaceI]];
label nFacePoints = f.size();
label nQuads = 0;
label nTris = 0;
f.nTrianglesQuads(mesh.points(), nTris, nQuads);
label nQuads = (nFacePoints - 2)/2;
label nTris = (nFacePoints - 2)%2;
nAddCells += nQuads + nTris;
}
@ -201,8 +200,8 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
}
else if (cellModel == prism)
{
// VTK has a different node order - their triangles point outwards!
// VTK has a different node order for VTK_WEDGE
// their triangles point outwards!
nodeIds[0] = cellShape[0];
nodeIds[1] = cellShape[2];
nodeIds[2] = cellShape[1];
@ -349,29 +348,34 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
// Mapping from additional point to cell
addPointCellLabels[addPointI] = cellI;
// Insert the new vertex from the cell-centre
label newVertexLabel = mesh.nPoints() + addPointI;
// The new vertex from the cell-centre
const label newVertexLabel = mesh.nPoints() + addPointI;
vtkInsertNextOpenFOAMPoint(vtkpoints, mesh.C()[cellI]);
// Whether to insert cell in place of original or not.
bool substituteCell = true;
const labelList& cFaces = mesh.cells()[cellI];
forAll(cFaces, cFaceI)
{
const face& f = mesh.faces()[cFaces[cFaceI]];
const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
label nFacePoints = f.size();
// Number of triangles and quads in decomposition
label nTris = 0;
label nQuads = 0;
f.nTrianglesQuads(mesh.points(), nTris, nQuads);
label nQuads = (nFacePoints - 2)/2;
label nTris = (nFacePoints - 2)%2;
// Do actual decomposition into triFcs and quadFcs.
faceList triFcs(nTris);
faceList quadFcs(nQuads);
label trii = 0;
label quadi = 0;
f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs);
label qpi = 0;
for (label quadi=0; quadi<nQuads; quadi++)
forAll(quadFcs, quadI)
{
label thisCellI = -1;
label thisCellI;
if (substituteCell)
{
@ -384,10 +388,29 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
superCells[addCellI++] = cellI;
}
nodeIds[0] = f[0];
nodeIds[1] = f[qpi + 1];
nodeIds[2] = f[qpi + 2];
nodeIds[3] = f[qpi + 3];
const face& quad = quadFcs[quadI];
// Ensure we have the correct orientation for the
// base of the primitive cell shape.
// If the cell is face owner, the orientation needs to be
// flipped.
// At the moment, VTK doesn't actually seem to care if
// negative cells are defined, but we'll do it anyhow
// (for safety).
if (isOwner)
{
nodeIds[0] = quad[3];
nodeIds[1] = quad[2];
nodeIds[2] = quad[1];
nodeIds[3] = quad[0];
}
else
{
nodeIds[0] = quad[0];
nodeIds[1] = quad[1];
nodeIds[2] = quad[2];
nodeIds[3] = quad[3];
}
nodeIds[4] = newVertexLabel;
vtkmesh->InsertNextCell
(
@ -395,13 +418,11 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
5,
nodeIds
);
qpi += 2;
}
if (nTris)
forAll(triFcs, triI)
{
label thisCellI = -1;
label thisCellI;
if (substituteCell)
{
@ -414,10 +435,23 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
superCells[addCellI++] = cellI;
}
nodeIds[0] = f[0];
nodeIds[1] = f[qpi + 1];
nodeIds[2] = f[qpi + 2];
const face& tri = triFcs[triI];
// See note above about the orientation.
if (isOwner)
{
nodeIds[0] = tri[2];
nodeIds[1] = tri[1];
nodeIds[2] = tri[0];
}
else
{
nodeIds[0] = tri[0];
nodeIds[1] = tri[1];
nodeIds[2] = tri[2];
}
nodeIds[3] = newVertexLabel;
vtkmesh->InsertNextCell
(
VTK_TETRA,