From a0d79333608624885837d818161794d4f3b2fd0c Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Fri, 17 Sep 2021 18:46:30 +0200 Subject: [PATCH] ENH: support vtu mesh subsetting, creation from cellShapes - support simple mesh subsetting to vtu formats to enable debug output for a subsection of an existing polyMesh - rudimentary support for VTK from cellShapes is intended for handling basic primitive cell shapes without a full blown polyMesh description. For example, // Create an empty polyMesh with points only polyMesh debugMesh ( io, std::move(points), faceList(), // faces labelList(), // owner labelList(), // neighbour true // syncPar ); // Establish the appropriate VTK sizing for the TET shapes vtk::vtuCells vtuCells; vtuCells.resetShapes(debugCutTets); vtuCells.nPoints(debugMesh.nPoints()); NOTE Since the vtk::internalMeshWriter only uses the polyMesh reference for the points, it is also possible to create the vtuCells description without a pointField (or from a different mesh description) and write out the connectivity using the pointField from a different mesh. --- .../test/foamCellZoneToVTK/Make/files | 3 + .../test/foamCellZoneToVTK/Make/options | 7 + .../foamCellZoneToVTK/foamCellZoneToVTK.C | 161 ++++ .../test/foamMeshToTet-vtk/Make/files | 4 + .../test/foamMeshToTet-vtk/Make/options | 7 + .../foamMeshToTet-vtk/foamMeshToTet-vtk.C | 96 +++ .../test/foamMeshToTet-vtk/writeVTKtetMesh.C | 204 +++++ src/fileFormats/vtk/part/foamVtkMeshMaps.H | 25 +- src/fileFormats/vtk/part/foamVtkMeshMapsI.H | 32 +- src/fileFormats/vtk/part/foamVtuCells.C | 193 ++++- src/fileFormats/vtk/part/foamVtuCells.H | 49 +- src/fileFormats/vtk/part/foamVtuCellsI.H | 20 +- src/fileFormats/vtk/part/foamVtuSizing.C | 355 ++++++++- src/fileFormats/vtk/part/foamVtuSizing.H | 147 +++- src/fileFormats/vtk/part/foamVtuSizingI.H | 37 +- ...uSizingTemplates.C => foamVtuSizingImpl.C} | 718 ++++++++++-------- 16 files changed, 1599 insertions(+), 459 deletions(-) create mode 100644 applications/test/foamCellZoneToVTK/Make/files create mode 100644 applications/test/foamCellZoneToVTK/Make/options create mode 100644 applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C create mode 100644 applications/test/foamMeshToTet-vtk/Make/files create mode 100644 applications/test/foamMeshToTet-vtk/Make/options create mode 100644 applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C create mode 100644 applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C rename src/fileFormats/vtk/part/{foamVtuSizingTemplates.C => foamVtuSizingImpl.C} (68%) diff --git a/applications/test/foamCellZoneToVTK/Make/files b/applications/test/foamCellZoneToVTK/Make/files new file mode 100644 index 0000000000..a4be1f0e87 --- /dev/null +++ b/applications/test/foamCellZoneToVTK/Make/files @@ -0,0 +1,3 @@ +foamCellZoneToVTK.C + +EXE = $(FOAM_USER_APPBIN)/foamCellZoneToVTK diff --git a/applications/test/foamCellZoneToVTK/Make/options b/applications/test/foamCellZoneToVTK/Make/options new file mode 100644 index 0000000000..d820876f44 --- /dev/null +++ b/applications/test/foamCellZoneToVTK/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -I$(LIB_SRC)/fileFormats/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfileFormats \ + -lmeshTools diff --git a/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C b/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C new file mode 100644 index 0000000000..a2f265db03 --- /dev/null +++ b/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C @@ -0,0 +1,161 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Application + foamCellZoneToVTK.C + +Description + Write tet-decomposed OpenFOAM mesh in VTK format. + For diagnostic purposes. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "timeSelector.H" +#include "Time.H" +#include "polyMesh.H" +#include "foamVtkInternalMeshWriter.H" + +using namespace Foam; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::addNote + ( + "Write OpenFOAM cellZone mesh to VTK" + ); + argList::addOption + ( + "cellZone", + "name", + "Convert mesh subset corresponding to specified cellZone" + ); + argList::addBoolOption + ( + "list", + "List names of cellZones and exit" + ); + timeSelector::addOptions(); + + #include "setRootCase.H" + + word cellZoneName; + args.readIfPresent("cellZone", cellZoneName); + + const bool optList = args.found("list"); + + if (optList) + { + if (!cellZoneName.empty()) + { + Info<< "specify -list or -cellZone, but not both!" << nl; + return 1; + } + } + else if (cellZoneName.empty()) + { + Info<< "Did not specify a cellZone!!" << nl; + return 1; + } + + #include "createTime.H" + + instantList timeDirs = timeSelector::select0(runTime, args); + + fileName exportName("zonemesh-" + cellZoneName); + + #include "createPolyMesh.H" + + forAll(timeDirs, timei) + { + runTime.setTime(timeDirs[timei], timei); + + polyMesh::readUpdateState state = mesh.readUpdate(); + + if (!timei || state != polyMesh::UNCHANGED) + { + fileName meshName(exportName); + if (state != polyMesh::UNCHANGED) + { + meshName += '_' + runTime.timeName(); + } + + if (optList) + { + Info<< "cellZones:" << nl; + for (const word& name : mesh.cellZones().names()) + { + Info<< " " << name << nl; + } + } + else + { + const cellZone* zonePtr = + mesh.cellZones().cfindZone(cellZoneName); + + Info<< "cellZone " << cellZoneName; + if (!zonePtr) + { + Info<< " ... not found" << nl; + continue; + } + Info<< nl; + + const cellZone& zn = *zonePtr; + + // Define a subset + vtk::vtuCells vtuCells; + vtuCells.reset(mesh, zn); + + vtk::internalMeshWriter writer + ( + mesh, + vtuCells, + fileName + ( + mesh.time().globalPath() / meshName + ) + ); + + writer.writeGeometry(); + + writer.beginCellData(); + writer.writeProcIDs(); + + Info<< "Wrote " << writer.output().name() << nl; + } + } + } + + Info<< "\nEnd\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/foamMeshToTet-vtk/Make/files b/applications/test/foamMeshToTet-vtk/Make/files new file mode 100644 index 0000000000..576d08e8c1 --- /dev/null +++ b/applications/test/foamMeshToTet-vtk/Make/files @@ -0,0 +1,4 @@ +foamMeshToTet-vtk.C +writeVTKtetMesh.C + +EXE = $(FOAM_USER_APPBIN)/foamMeshToTet-vtk diff --git a/applications/test/foamMeshToTet-vtk/Make/options b/applications/test/foamMeshToTet-vtk/Make/options new file mode 100644 index 0000000000..d820876f44 --- /dev/null +++ b/applications/test/foamMeshToTet-vtk/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -I$(LIB_SRC)/fileFormats/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfileFormats \ + -lmeshTools diff --git a/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C b/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C new file mode 100644 index 0000000000..d1d078949e --- /dev/null +++ b/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C @@ -0,0 +1,96 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Application + foamMeshToTet-vtk + +Description + Write tet-decomposed OpenFOAM mesh in VTK format. + For diagnostic purposes. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "timeSelector.H" +#include "Time.H" +#include "polyMesh.H" + +namespace Foam +{ + void writeVTKtetMesh(const fileName& output, const polyMesh& mesh); +} + +using namespace Foam; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::addNote + ( + "Write tet-decomposed OpenFOAM mesh in VTK" + ); + argList::noParallel(); + timeSelector::addOptions(); + + #include "setRootCase.H" + #include "createTime.H" + + instantList timeDirs = timeSelector::select0(runTime, args); + + fileName exportName = "tetmesh"; + if (args.found("case")) + { + exportName += '-' + args.globalCaseName(); + } + + #include "createPolyMesh.H" + + forAll(timeDirs, timei) + { + runTime.setTime(timeDirs[timei], timei); + + polyMesh::readUpdateState state = mesh.readUpdate(); + + if (!timei || state != polyMesh::UNCHANGED) + { + fileName meshName(exportName); + if (state != polyMesh::UNCHANGED) + { + meshName += '_' + runTime.timeName(); + } + + writeVTKtetMesh(meshName, mesh); + } + } + + Info<< "\nEnd\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C b/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C new file mode 100644 index 0000000000..ea7157fac0 --- /dev/null +++ b/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C @@ -0,0 +1,204 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "polyMesh.H" +#include "Fstream.H" +#include "tetMatcher.H" +#include "foamVtkInternalMeshWriter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +void writeVTKtetMesh(const fileName& output, const polyMesh& mesh_) +{ + #if WM_LABEL_SIZE == 64 + # define FOAM_VTK_LABEL_NAME "vtktypeint64" + #else + # define FOAM_VTK_LABEL_NAME "int" + #endif + + const faceList& faces = mesh_.faces(); + const labelList& faceOwner = mesh_.faceOwner(); + + label nTets = 0; + + for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei) + { + nTets += 2 * faces[facei].nTriangles(); + } + for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei) + { + nTets += faces[facei].nTriangles(); + } + const label nPoints = (mesh_.nPoints() + mesh_.nCells()); + + OFstream os(output + ".vtk"); + + Info<< "Write: " << os.name() << nl; + + os << "# vtk DataFile Version 5.1" << nl + << "tet-mesh" << nl + << "ASCII" << nl + << "DATASET UNSTRUCTURED_GRID" << nl + << nl; + + os << "POINTS " << nPoints << " float" << nl; + for (const point& p : mesh_.points()) + { + os << float(p.x()) << ' ' + << float(p.y()) << ' ' + << float(p.z()) << nl; + } + for (const point& p : mesh_.cellCentres()) + { + os << float(p.x()) << ' ' + << float(p.y()) << ' ' + << float(p.z()) << nl; + } + os << nl; + + os << "CELLS " << (1 + nTets) << ' ' << (4 * nTets) << nl; + + os << "OFFSETS " << FOAM_VTK_LABEL_NAME << nl + << 0; // begin offset = 0 + { + label offset = 0; + for (label teti = 0; teti < nTets; ++teti) + { + offset += 4; + os << ' ' << offset; + } + os << nl << nl; + } + + labelList nLocalTets(mesh_.nCells(), Zero); + + os << nl + << "CONNECTIVITY " << FOAM_VTK_LABEL_NAME << nl; + + for (label celli = 0; celli < mesh_.nCells(); ++celli) + { + const cell& cFaces = mesh_.cells()[celli]; + + if (tetMatcher::test(mesh_, celli)) + { + // Tet: no cell-centre decomposition + + const label facei = cFaces[0]; + const face& f0 = faces[facei]; + + // Get the other point from f1. Tbd: check if not duplicate face + // (ACMI / ignoreBoundaryFaces_). + const face& f1 = faces[cFaces[1]]; + label apexi = -1; + forAll(f1, fp) + { + apexi = f1[fp]; + if (!f0.found(apexi)) + { + break; + } + } + + const label p0 = f0[0]; + label p1 = f0[1]; + label p2 = f0[2]; + + if (faceOwner[facei] == celli) + { + std::swap(p1, p2); + } + + ++nLocalTets[celli]; + os << p0 << ' ' << p1 << ' ' << p2 << ' ' << apexi << nl; + } + else + { + for (const label facei : cFaces) + { + const face& f = faces[facei]; + + label fp0 = mesh_.tetBasePtIs()[facei]; + + // Fallback + if (fp0 < 0) + { + fp0 = 0; + } + + const label p0 = f[fp0]; + label fp = f.fcIndex(fp0); + for (label i = 2; i < f.size(); ++i) + { + label p1 = f[fp]; + fp = f.fcIndex(fp); + label p2 = f[fp]; + + if (faceOwner[facei] == celli) + { + std::swap(p1, p2); + } + + ++nLocalTets[celli]; + os << p0 << ' ' << p1 << ' ' << p2 << ' ' + << (mesh_.nPoints() + celli) << nl; + } + } + } + } + + os << nl + << "CELL_TYPES " << nTets << nl; + + for (label teti = 0; teti < nTets; ++teti) + { + if (teti) os << ' '; + os << vtk::cellType::VTK_TETRA; + } + os << nl; + + os << nl << "CELL_DATA " << nTets << nl + << "FIELD FieldData " << 1 << nl; + + os << "cellID " << 1 << ' ' << nTets << " int" << nl; + forAll(nLocalTets, celli) + { + label n = nLocalTets[celli]; + while (n--) + { + os << ' ' << celli; + } + os << nl; + } +} + +} // End namespace Foam + + +// ************************************************************************* // diff --git a/src/fileFormats/vtk/part/foamVtkMeshMaps.H b/src/fileFormats/vtk/part/foamVtkMeshMaps.H index c2c54cc049..474c2d158b 100644 --- a/src/fileFormats/vtk/part/foamVtkMeshMaps.H +++ b/src/fileFormats/vtk/part/foamVtkMeshMaps.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017-2020 OpenCFD Ltd. + Copyright (C) 2017-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -75,11 +75,8 @@ public: // Constructors - //- Default construct: zero-sized, no reserved size - inline foamVtkMeshMaps(); - - //- Construct with reserved size - inline explicit foamVtkMeshMaps(const label size); + //- Default construct: zero-sized, no reserved sizes + foamVtkMeshMaps() = default; // Member Functions @@ -87,26 +84,26 @@ public: // Access //- Original cell ids for all cells (regular and decomposed). - // A regular mesh comprising only primitive cell types, this will just + // For regular mesh comprising only primitive cell types, this will // be an identity list. However, for subsetted meshes and decomposed // cells this becomes a useful means of mapping from the original mesh. - inline const labelList& cellMap() const; + inline const labelList& cellMap() const noexcept; //- Write access to original cell ids - inline DynamicList