65#define SNAP_END_ENCODE(pos, val) (((val) << (4 + 2 * pos)))
66#define SNAP_END_VALUE(pos, val) (((val) >> (4 + 2 * pos)) & 0x3)
94 if (doSnap && cutIndex && cutIndex != 0xF)
99 p0 -= val;
if (cutIndex & 1)
p0 *= -1;
100 p1 -= val;
if (cutIndex & 2) p1 *= -1;
101 p2 -= val;
if (cutIndex & 4) p2 *= -1;
102 p3 -= val;
if (cutIndex & 8) p3 *= -1;
107 #undef ADD_SNAP_INDEX
108 #define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2) \
109 switch (cutIndex & (idx1 | idx2)) \
113 cutIndex |= SNAP_END_ENCODE \
116 ((d1 * 100 < d2) ? 1 : (d2 * 100 < d1) ? 2 : 0) \
127 #undef ADD_SNAP_INDEX
145 if (a !=
b &&
b != c && c != a)
182Foam::isoSurfaceTopo::tetCutAddressing::tetCutAddressing
184 const label nCutCells,
186 const bool useDebugCuts
189 vertsToPointLookup_(12*nCutCells),
192 pointToFace_(10*nCutCells),
193 pointFromDiag_(10*nCutCells),
195 pointToVerts_(10*nCutCells),
196 cutPoints_(12*nCutCells),
199 debugCutTetsOn_(useDebugCuts)
210 snapVertsLookup_.resize(4*nCutCells);
214 debugCutTets_.reserve(6*nCutCells);
221void Foam::isoSurfaceTopo::tetCutAddressing::clearDebug()
223 debugCutTets_.clearStorage();
227void Foam::isoSurfaceTopo::tetCutAddressing::clearDiagonal()
229 pointToFace_.clearStorage();
230 pointFromDiag_.clearStorage();
234void Foam::isoSurfaceTopo::tetCutAddressing::clearHashes()
236 vertsToPointLookup_.clear();
237 snapVertsLookup_.clear();
241Foam::label Foam::isoSurfaceTopo::tetCutAddressing::generatePoint
252 label pointi = vertsToPointLookup_.lookup(vertices, -1);
255 bool addNewPoint(
true);
257 const label snapPointi =
260 : (snapEnd == 2) ?
vertices.second()
264 if (snapPointi == -1)
267 pointi = pointToVerts_.size();
268 pointToVerts_.append(vertices);
274 edgeIsDiagonal =
false;
276 pointi = snapVertsLookup_.lookup(snapPointi, -1);
277 addNewPoint = (pointi == -1);
280 pointi = pointToVerts_.size();
281 snapVertsLookup_.insert(snapPointi, pointi);
282 pointToVerts_.append(edge(snapPointi, snapPointi));
288 pointToFace_.append(facei);
289 pointFromDiag_.append(edgeIsDiagonal);
292 vertsToPointLookup_.insert(vertices, pointi);
299bool Foam::isoSurfaceTopo::tetCutAddressing::generatePoints
302 const int tetCutIndex,
303 const tetCell& tetLabels,
306 const FixedList<bool, 6>& edgeIsDiagonal
310 const label nCutPointsOld(cutPoints_.size());
313 switch (tetCutIndex & 0x0F)
320 case 0x0E: flip =
true; [[fallthrough]];
359 case 0x0D: flip =
true; [[fallthrough]];
398 case 0x0C: flip =
true; [[fallthrough]];
448 case 0x0B: flip =
true; [[fallthrough]];
487 case 0x0A: flip =
true; [[fallthrough]];
537 case 0x09: flip =
true; [[fallthrough]];
587 case 0x07: flip =
true; [[fallthrough]];
626 const bool added(nCutPointsOld != cutPoints_.size());
628 if (added && debugCutTetsOn_)
630 debugCutTets_.append(tetLabels.shape());
640void Foam::isoSurfaceTopo::generateTriPoints
644 const labelList& tetBasePtIs,
645 tetCutAddressing& tetCutAddr
651 const bool doSnap = this->
snap();
658 const label facei = cFaces[0];
659 const face& f0 = faces[facei];
663 const face& f1 = faces[cFaces[1]];
668 if (!f0.found(apexi))
674 const label
p0 = f0[0];
678 if (faceOwner[facei] == celli)
683 const tetCell tetLabels(
p0, p1, p2, apexi);
684 const int tetCutIndex
697 tetCutAddr.generatePoints
702 FixedList<bool, 6>(
false)
707 for (
const label facei : cFaces)
718 const face&
f = faces[facei];
720 label fp0 = tetBasePtIs[facei];
728 const label
p0 =
f[fp0];
730 for (label i = 2; i <
f.
size(); ++i)
736 FixedList<bool, 6> edgeIsDiagonal(
false);
737 if (faceOwner[facei] == celli)
740 if (i != 2) edgeIsDiagonal[1] =
true;
741 if (i !=
f.
size()-1) edgeIsDiagonal[0] =
true;
745 if (i != 2) edgeIsDiagonal[0] =
true;
746 if (i !=
f.
size()-1) edgeIsDiagonal[1] =
true;
750 const int tetCutIndex
763 tetCutAddr.generatePoints
778void Foam::isoSurfaceTopo::triangulateOutside
780 const bool filterDiag,
787 DynamicList<face>& compactFaces,
788 DynamicList<label>& compactCellIDs
805 compactFaces.append(face(loop.size()));
806 face&
f = compactFaces.
last();
811 const label pointi =
mp[loop[i]];
812 if (filterDiag && pointFromDiag[pointi])
814 const label prevPointi =
mp[loop[loop.fcIndex(i)]];
817 pointFromDiag[prevPointi]
818 && (pointToFace[pointi] != pointToFace[prevPointi])
843 const label pointi =
mp[loop[i]];
847 compactCellIDs.
append(cellID);
853void Foam::isoSurfaceTopo::removeInsidePoints
856 const bool filterDiag,
864 DynamicList<label>& compactCellIDs
869 compactCellIDs.clear();
870 compactCellIDs.reserve(
s.size()/4);
872 DynamicList<face> compactFaces(
s.size()/4);
874 for (label celli = 0; celli < start.size()-1; ++celli)
878 const label nTris = start[celli+1]-start[celli];
884 SubList<face>(
s, nTris, start[celli]),
902 s.swapFaces(compactFaces);
928 Pout<<
"isoSurfaceTopo:" <<
nl
931 <<
" isoValue : " << iso <<
nl
935 <<
" ignoreCells : " << ignoreCells.
count()
942 label nBlockedCells = 0;
945 nBlockedCells +=
blockCells(cellCutType_, ignoreCells);
963 Pout<<
"isoSurfaceTopo : candidate cells cut "
965 <<
" blocked " << nBlockedCells
969 if (debug && isA<fvMesh>(
mesh))
971 const auto& fvmesh = dynamicCast<const fvMesh>(
mesh);
977 "isoSurfaceTopo.cutType",
978 fvmesh.time().timeName(),
990 forAll(cellCutType_, celli)
992 debugFld[celli] = cellCutType_[celli];
1008 forAll(cellCutType_, celli)
1010 if ((cellCutType_[celli] & realCut) != 0)
1012 debugCutCells[nout] = celli;
1014 if (nout >= nCutCells)
break;
1029 / (
"isoSurfaceTopo." + timeDesc +
"-cutCells")
1044 Info<<
"isoSurfaceTopo : (debug) wrote "
1052 tetCutAddressing tetCutAddr
1060 for (label celli = 0; celli <
mesh_.
nCells(); ++celli)
1062 startTri[celli] = tetCutAddr.nFaces();
1076 startTri.
last() = tetCutAddr.nFaces();
1079 tetBasePtIs.
clear();
1080 tetCutAddr.clearHashes();
1086 auto& verts = tetCutAddr.cutPoints();
1089 for (
face&
f : allTriFaces)
1092 f[0] = verts[verti++];
1093 f[1] = verts[verti++];
1094 f[2] = verts[verti++];
1096 verts.clearStorage();
1102 for (label celli = 0; celli < startTri.
size()-1; ++celli)
1108 (startTri[celli+1] - startTri[celli]),
1114 pointToVerts_.
transfer(tetCutAddr.pointToVerts());
1129 std::move(allTriPoints),
1130 std::move(allTriFaces),
1136 Pout<<
"isoSurfaceTopo : generated "
1142 if ((debug & 8) && (params.
filter() != filterType::NONE))
1144 const Mesh&
s = *
this;
1153 / (
"isoSurfaceTopo." + timeDesc +
"-triangles")
1174 const edge& verts = pointToVerts_[i];
1180 if (tetCutAddr.pointFromDiag().test(i))
1183 pointStatus[i] = -1;
1187 writer.write(
"point-status", pointStatus);
1190 Info<<
"isoSurfaceTopo : (debug) wrote "
1205 if (params.
filter() == filterType::NONE)
1213 s.compactPoints(pointMap);
1232 params.
filter() == filterType::DIAGCELL
1233 || params.
filter() == filterType::NONMANIFOLD
1235 tetCutAddr.pointFromDiag(),
1236 tetCutAddr.pointToFace(),
1242 s.compactPoints(pointMap);
1249 Pout<<
"isoSurfaceTopo :"
1250 " after removing cell centre and face-diag triangles: "
1258 tetCutAddr.clearDiagonal();
1267 (params.
filter() == filterType::NONMANIFOLD)
1268 || tetCutAddr.debugCutTetsOn()
1309 cellCutType_.
clear();
1313 if (tetCutAddr.debugCutTetsOn())
1317 const auto& debugCuts = tetCutAddr.debugCutTets();
1334 / (
"isoSurfaceTopo." + timeDesc +
"-cutTets")
1347 forAll(cutTetQuality, teti)
1357 writer.writeCellData(
"tetQuality", cutTetQuality);
1367 for (
const edge& verts : pointToVerts_)
1369 if (verts.first() == verts.second())
1372 pointStatus[verts.first()] = 1;
1376 writer.writePointData(
"point-status", pointStatus);
1379 Info<<
"isoSurfaceTopo : (debug) wrote "
1390 const Mesh&
s = *
this;
1395 openEdgeIds.
resize(2*
s.size());
1400 if (eFaces.
size() == 1)
1404 const edge&
e = surfEdges[edgei];
1405 const edge& verts0 = pointToVerts_[mp[
e.
first()]];
1406 const edge& verts1 = pointToVerts_[mp[
e.second()]];
1412 && isProtectedPoint.
test(verts1.
first())
1421 openEdgeIds.
insert(edgei);
1426 const label nOpenEdges
1446 / (
"isoSurfaceTopo." + timeDesc +
"-openEdges")
1456 Info<<
"isoSurfaceTopo : (debug) wrote "
1457 << nOpenEdges <<
" open edges: "
1462 Info<<
"isoSurfaceTopo : no open edges" <<
nl;
1469 const Mesh&
s = *
this;
1478 / (
"isoSurfaceTopo." + timeDesc +
"-surface")
1499 const edge& verts = pointToVerts_[i];
1507 writer.write(
"point-status", pointStatus);
1510 Info<<
"isoSurfaceTopo : (debug) wrote "
1516 tetCutAddr.clearDebug();
1519 if (params.
filter() == filterType::NONMANIFOLD)
1551 label nEdgeRemove = 0;
1556 if (eFaces.
size() == 1)
1560 const edge&
e = surfEdges[edgei];
1561 const edge& verts0 = pointToVerts_[mp[
e.
first()]];
1562 const edge& verts1 = pointToVerts_[mp[
e.second()]];
1568 && isProtectedPoint.
test(verts1.
first())
1585 Pout<<
"isoSurfaceTopo :"
1588 <<
" faces on " << nEdgeRemove <<
" open edges" <<
endl;
1604 if (surf.
size() != faceAddr.
size())
1624 Mesh::subsetMesh(include, pointMap,
faceMap)
1626 Mesh::transfer(filtered);
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
T & first() noexcept
The first element of the list, position [0].
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
void resize(const label sz)
Resize the hash table for efficiency.
label size() const noexcept
The number of elements in table.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
fileName objectRelPath() const
The object path relative to the root.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void transfer(List< T > &list)
SubList< label > subList
Declare type of subList.
void append(const T &val)
Append an element at the end of the list.
void resize(const label len)
Adjust allocated size of list.
void clear()
Clear the list, i.e. set size to zero.
label size() const
The surface size is the number of faces.
void resize(const label numElem, const unsigned int val=0u)
Reset addressable list size, does not shrink the allocated size.
const T & second() const noexcept
Return second element, which is also the last element.
A list of faces which address into the list of points.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
fileName globalPath() const
Return global path for the case.
label timeIndex() const noexcept
Return current time index.
A List with indirect addressing. Like IndirectList but does not store addressing.
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
label fcIndex(const label i) const noexcept
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
unsigned int count(const bool on=true) const
Count number of bits set.
void set(const bitSet &bitset)
Set specified bits from another bitset.
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
scalar mag() const
The magnitude of the bounding box span.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Face selection method for createBaffles.
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
Low-level components common to various iso-surface algorithms.
bitSet ignoreBoundaryFaces_
Optional boundary faces to ignore.
@ BLOCKED
Blocked (never cut)
@ ANYCUT
Any cut type (bitmask)
@ TETCUT
Cell cut is a tet.
const scalarField & pVals_
Point values.
label calcCellCuts(List< cutType > &cuts) const
Populate a list of candidate cell cuts using getCellCutType()
labelList meshCells_
For every face, the original cell in mesh.
const scalar iso_
Iso value.
void ignoreCyclics()
Set ignoreBoundaryFaces to ignore cyclics (cyclicACMI)
const polyMesh & mesh() const noexcept
The mesh for which the iso-surface is associated.
const polyMesh & mesh_
Reference to mesh.
label blockCells(UList< cutType > &cuts, const bitSet &ignoreCells) const
Mark ignoreCells as BLOCKED.
const scalarField & cVals_
Cell values.
Preferences for controlling iso-surface algorithms.
filterType filter() const noexcept
Get current filter type.
static const Enum< filterType > filterNames
Names for the filtering types.
const boundBox & getClipBounds() const noexcept
Get optional clipping bounding box.
bool snap() const noexcept
Get point snapping flag.
Marching tet iso surface algorithm with optional filtering to keep only points originating from mesh ...
void inplaceSubsetMesh(const bitSet &include)
Subset the surface using the selected faces.
tmp< Field< Type > > interpolateTemplate(const Field< Type > &cellData, const Field< Type > &pointData) const
Interpolates cellData and pointData fields.
const Time & time() const noexcept
Return time registry.
static labelList adjustTetBasePtIs(const polyMesh &mesh, const bool report=false)
Return an adjusted list of tet base points.
Mesh consisting of general polyhedral cells.
virtual const faceList & faces() const
Return raw faces.
const boundBox & bounds() const
Return mesh bounding box.
virtual const labelList & faceOwner() const
Return face owner.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
virtual bool write(const bool valid=true) const
Write using setting from DB.
@ OUTSIDE
A location outside the volume.
Write an OpenFOAM volume (internal) geometry and internal fields as a vtu file or a legacy vtk file.
virtual bool beginPointData(label nFields=0)
Begin PointData for specified number of fields.
Write edge/points (optionally with fields) as a vtp file or a legacy vtk file.
Write faces/points (optionally with fields) as a vtp file or a legacy vtk file.
A deep-copy description of an OpenFOAM volume mesh in data structures suitable for VTK UnstructuredGr...
void resetShapes(const UList< cellShape > &shapes)
Reset sizing using primitive shapes only (ADVANCED USAGE)
void addPointCellLabels(const labelUList &cellIds)
Define which additional cell-centres are to be used (ADVANCED!)
void reset(const polyMesh &mesh)
void setNumPoints(label n) noexcept
Alter number of mesh points (ADVANCED USAGE)
label nFieldPoints() const noexcept
Number of field points = nPoints + nAddPoints.
A class for handling words, derived from Foam::string.
static word printf(const char *fmt, const PrimitiveType &val)
Use a printf-style formatter for a primitive.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & p0
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Convenience macros for instantiating iso-surface interpolate methods.
#define defineIsoSurfaceInterpolateMethods(ThisClass)
#define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2)
#define SNAP_END_VALUE(pos, val)
const dimensionedScalar mp
Proton mass.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
const dimensionSet dimless
Dimensionless.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
List< label > labelList
A List of labels.
static void appendTriLabels(DynamicList< label > &verts, const label a, const label b, const label c, const bool flip)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
MinMax< label > minMax(const labelHashSet &set)
Find the min/max values of labelHashSet.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< labelList > labelListList
A List of labelList.
pointField vertices(const blockVertexList &bvl)
static constexpr const zero Zero
Global zero (0)
UList< bool > boolUList
A UList of bools.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< surfZone > surfZoneList
List< face > faceList
A List of faces.
static const point & getMeshPointRef(const polyMesh &mesh, const label pointi)
UList< label > labelUList
A UList of labels.
static int getTetCutIndex(scalar p0, scalar p1, scalar p2, scalar p3, const scalar val, const bool doSnap) noexcept
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.