63 #undef SNAP_END_ENCODE 64 #define SNAP_END_ENCODE(pos, val) (((val) << (4 + 2 * pos))) 65 #define SNAP_END_VALUE(pos, val) (((val) >> (4 + 2 * pos)) & 0x3) 93 if (doSnap && cutIndex && cutIndex != 0xF)
98 p0 -= val;
if (cutIndex & 1)
p0 *= -1;
99 p1 -= val;
if (cutIndex & 2) p1 *= -1;
100 p2 -= val;
if (cutIndex & 4) p2 *= -1;
101 p3 -= val;
if (cutIndex & 8) p3 *= -1;
106 #undef ADD_SNAP_INDEX 107 #define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2) \ 108 switch (cutIndex & (idx1 | idx2)) \ 112 cutIndex |= SNAP_END_ENCODE \ 115 ((d1 * 100 < d2) ? 1 : (d2 * 100 < d1) ? 2 : 0) \ 126 #undef ADD_SNAP_INDEX 137 DynamicList<label>& verts,
144 if (a !=
b &&
b !=
c &&
c != a)
164 const polyMesh&
mesh,
181 Foam::isoSurfaceTopo::tetCutAddressing::tetCutAddressing
183 const label nCutCells,
185 const bool useDebugCuts
188 vertsToPointLookup_(12*nCutCells),
191 pointToFace_(10*nCutCells),
192 pointFromDiag_(10*nCutCells),
194 pointToVerts_(10*nCutCells),
195 cutPoints_(12*nCutCells),
198 debugCutTetsOn_(useDebugCuts)
209 snapVertsLookup_.reserve(2*nCutCells);
213 debugCutTets_.reserve(6*nCutCells);
220 void Foam::isoSurfaceTopo::tetCutAddressing::clearDebug()
222 debugCutTets_.clearStorage();
226 void Foam::isoSurfaceTopo::tetCutAddressing::clearDiagonal()
228 pointToFace_.clearStorage();
229 pointFromDiag_.clearStorage();
233 void Foam::isoSurfaceTopo::tetCutAddressing::clearHashes()
235 vertsToPointLookup_.clear();
236 snapVertsLookup_.clear();
240 Foam::label Foam::isoSurfaceTopo::tetCutAddressing::generatePoint
251 label pointi = vertsToPointLookup_.lookup(
vertices, -1);
254 bool addNewPoint(
true);
256 const label snapPointi =
259 : (snapEnd == 2) ?
vertices.second()
263 if (snapPointi == -1)
266 pointi = pointToVerts_.size();
273 edgeIsDiagonal =
false;
275 pointi = snapVertsLookup_.lookup(snapPointi, -1);
276 addNewPoint = (pointi == -1);
279 pointi = pointToVerts_.size();
280 snapVertsLookup_.insert(snapPointi, pointi);
281 pointToVerts_.append(edge(snapPointi, snapPointi));
287 pointToFace_.append(facei);
288 pointFromDiag_.append(edgeIsDiagonal);
291 vertsToPointLookup_.insert(
vertices, pointi);
298 bool Foam::isoSurfaceTopo::tetCutAddressing::generatePoints
301 const int tetCutIndex,
302 const tetCell& tetLabels,
305 const FixedList<bool, 6>& edgeIsDiagonal
309 const label nCutPointsOld(cutPoints_.size());
312 switch (tetCutIndex & 0x0F)
319 case 0x0E: flip =
true; [[fallthrough]];
358 case 0x0D: flip =
true; [[fallthrough]];
397 case 0x0C: flip =
true; [[fallthrough]];
447 case 0x0B: flip =
true; [[fallthrough]];
486 case 0x0A: flip =
true; [[fallthrough]];
536 case 0x09: flip =
true; [[fallthrough]];
586 case 0x07: flip =
true; [[fallthrough]];
625 const bool added(nCutPointsOld != cutPoints_.size());
627 if (added && debugCutTetsOn_)
629 debugCutTets_.append(tetLabels.shape());
639 void Foam::isoSurfaceTopo::generateTriPoints
644 tetCutAddressing& tetCutAddr
650 const bool doSnap = this->
snap();
657 const label facei = cFaces[0];
658 const face& f0 = faces[facei];
662 const face& f1 = faces[cFaces[1]];
667 if (!f0.found(apexi))
673 const label
p0 = f0[0];
677 if (faceOwner[facei] == celli)
682 const tetCell tetLabels(
p0, p1, p2, apexi);
683 const int tetCutIndex
696 tetCutAddr.generatePoints
701 FixedList<bool, 6>(
false)
706 for (
const label facei : cFaces)
717 const face&
f = faces[facei];
719 label fp0 = tetBasePtIs[facei];
727 const label
p0 =
f[fp0];
728 label fp =
f.fcIndex(fp0);
729 for (label i = 2; i <
f.size(); ++i)
735 FixedList<bool, 6> edgeIsDiagonal(
false);
736 if (faceOwner[facei] == celli)
739 if (i != 2) edgeIsDiagonal[1] =
true;
740 if (i !=
f.size()-1) edgeIsDiagonal[0] =
true;
744 if (i != 2) edgeIsDiagonal[0] =
true;
745 if (i !=
f.size()-1) edgeIsDiagonal[1] =
true;
749 const int tetCutIndex
762 tetCutAddr.generatePoints
777 void Foam::isoSurfaceTopo::triangulateOutside
779 const bool filterDiag,
786 DynamicList<face>& compactFaces,
787 DynamicList<label>& compactCellIDs
809 const label pointi =
mp[loop[i]];
810 if (filterDiag && pointFromDiag[pointi])
812 const label prevPointi =
mp[loop[loop.fcIndex(i)]];
815 pointFromDiag[prevPointi]
816 && (pointToFace[pointi] != pointToFace[prevPointi])
841 const label pointi =
mp[loop[i]];
845 compactCellIDs.
append(cellID);
851 void Foam::isoSurfaceTopo::removeInsidePoints
854 const bool filterDiag,
862 DynamicList<label>& compactCellIDs
867 compactCellIDs.clear();
868 compactCellIDs.reserve(
s.size()/4);
870 DynamicList<face> compactFaces(
s.size()/4);
872 for (label celli = 0; celli < start.size()-1; ++celli)
876 const label nTris = start[celli+1]-start[celli];
882 SubList<face>(
s, nTris, start[celli]),
900 s.swapFaces(compactFaces);
926 Pout<<
"isoSurfaceTopo:" <<
nl 929 <<
" isoValue : " << iso <<
nl 933 <<
" ignoreCells : " << ignoreCells.
count()
940 label nBlockedCells = 0;
943 nBlockedCells +=
blockCells(cellCutType_, ignoreCells);
961 Pout<<
"isoSurfaceTopo : candidate cells cut " 963 <<
" blocked " << nBlockedCells
969 const auto& fvmesh = refCast<const fvMesh>(
mesh);
975 "isoSurfaceTopo.cutType",
976 fvmesh.time().timeName(),
988 forAll(cellCutType_, celli)
990 debugFld[celli] = cellCutType_[celli];
993 Info<<
"Writing cut types: " << debugField.objectRelPath() <<
nl;
1002 constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
1006 forAll(cellCutType_, celli)
1008 if ((cellCutType_[celli] & realCut) != 0)
1010 debugCutCells[nout] = celli;
1012 if (nout >= nCutCells)
break;
1017 vtk::vtuCells vtuCells;
1018 vtuCells.reset(
mesh_, debugCutCells);
1020 vtk::internalMeshWriter
writer 1027 / (
"isoSurfaceTopo." + timeDesc +
"-cutCells")
1042 Info<<
"isoSurfaceTopo : (debug) wrote " 1050 tetCutAddressing tetCutAddr
1058 for (label celli = 0; celli <
mesh_.
nCells(); ++celli)
1060 startTri[celli] = tetCutAddr.nFaces();
1061 if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
1067 bool(cellCutType_[celli] & cutType::TETCUT),
1074 startTri.back() = tetCutAddr.nFaces();
1077 tetBasePtIs.clear();
1078 tetCutAddr.clearHashes();
1082 faceList allTriFaces(startTri.back());
1084 auto& verts = tetCutAddr.cutPoints();
1087 for (face&
f : allTriFaces)
1090 f[0] = verts[verti++];
1091 f[1] = verts[verti++];
1092 f[2] = verts[verti++];
1100 for (label celli = 0; celli < startTri.size()-1; ++celli)
1106 (startTri[celli+1] - startTri[celli]),
1112 pointToVerts_.
transfer(tetCutAddr.pointToVerts());
1127 std::move(allTriPoints),
1128 std::move(allTriFaces),
1134 Pout<<
"isoSurfaceTopo : generated " 1142 const Mesh&
s = *
this;
1144 vtk::surfaceWriter
writer 1151 / (
"isoSurfaceTopo." + timeDesc +
"-triangles")
1172 const edge& verts = pointToVerts_[i];
1173 if (verts.first() == verts.second())
1178 if (tetCutAddr.pointFromDiag().test(i))
1181 pointStatus[i] = -1;
1185 writer.write(
"point-status", pointStatus);
1188 Info<<
"isoSurfaceTopo : (debug) wrote " 1211 s.compactPoints(pointMap);
1212 pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
1223 DynamicList<label> compactCellIDs;
1233 tetCutAddr.pointFromDiag(),
1234 tetCutAddr.pointToFace(),
1240 s.compactPoints(pointMap);
1242 pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
1247 Pout<<
"isoSurfaceTopo :" 1248 " after removing cell centre and face-diag triangles: " 1256 tetCutAddr.clearDiagonal();
1262 bitSet isProtectedPoint;
1266 || tetCutAddr.debugCutTetsOn()
1307 cellCutType_.clear();
1311 if (tetCutAddr.debugCutTetsOn())
1315 const auto& debugCuts = tetCutAddr.debugCutTets();
1318 vtk::vtuCells vtuCells;
1319 vtuCells.resetShapes(debugCuts);
1325 vtk::internalMeshWriter
writer 1332 / (
"isoSurfaceTopo." + timeDesc +
"-cutTets")
1344 Field<scalar> cutTetQuality(debugCuts.size());
1345 forAll(cutTetQuality, teti)
1355 writer.writeCellData(
"tetQuality", cutTetQuality);
1365 for (
const edge& verts : pointToVerts_)
1367 if (verts.first() == verts.second())
1370 pointStatus[verts.first()] = 1;
1374 writer.writePointData(
"point-status", pointStatus);
1377 Info<<
"isoSurfaceTopo : (debug) wrote " 1388 const Mesh&
s = *
this;
1393 openEdgeIds.reserve(
s.size());
1398 if (eFaces.size() == 1)
1402 const edge&
e = surfEdges[edgei];
1403 const edge& verts0 = pointToVerts_[
mp[
e.first()]];
1404 const edge& verts1 = pointToVerts_[
mp[
e.second()]];
1408 isProtectedPoint.test(verts0.first())
1409 && isProtectedPoint.test(verts0.second())
1410 && isProtectedPoint.test(verts1.first())
1411 && isProtectedPoint.test(verts1.second())
1419 openEdgeIds.insert(edgei);
1424 const label nOpenEdges
1434 openEdgeIds.sortedToc()
1444 / (
"isoSurfaceTopo." + timeDesc +
"-openEdges")
1454 Info<<
"isoSurfaceTopo : (debug) wrote " 1455 << nOpenEdges <<
" open edges: " 1460 Info<<
"isoSurfaceTopo : no open edges" <<
nl;
1467 const Mesh& s = *
this;
1469 vtk::surfaceWriter
writer 1476 / (
"isoSurfaceTopo." + timeDesc +
"-surface")
1497 const edge& verts = pointToVerts_[i];
1498 if (verts.first() == verts.second())
1505 writer.write(
"point-status", pointStatus);
1508 Info<<
"isoSurfaceTopo : (debug) wrote " 1514 tetCutAddr.clearDebug();
1531 bitSet faceSelection;
1538 UIndirectList<face>(surf, faceAddr),
1542 faceSelection.clear();
1543 faceSelection.resize(erosion.size());
1546 const edgeList& surfEdges = erosion.edges();
1549 label nEdgeRemove = 0;
1554 if (eFaces.size() == 1)
1558 const edge&
e = surfEdges[edgei];
1559 const edge& verts0 = pointToVerts_[
mp[
e.first()]];
1560 const edge& verts1 = pointToVerts_[
mp[
e.second()]];
1564 isProtectedPoint.test(verts0.first())
1565 && isProtectedPoint.test(verts0.second())
1566 && isProtectedPoint.test(verts1.first())
1567 && isProtectedPoint.test(verts1.second())
1575 faceSelection.set(eFaces[0]);
1583 Pout<<
"isoSurfaceTopo :" 1584 <<
" removing " << faceSelection.count()
1585 <<
" / " << faceSelection.size()
1586 <<
" faces on " << nEdgeRemove <<
" open edges" <<
endl;
1602 if (surf.size() != faceAddr.size())
1604 faceSelection.clear();
1605 faceSelection.resize(surf.size());
1606 faceSelection.set(faceAddr);
1622 Mesh::subsetMesh(include, pointMap,
faceMap)
1624 Mesh::transfer(filtered);
1626 meshCells_ = UIndirectList<label>(meshCells_,
faceMap)();
1628 pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
void ignoreCyclics()
Set ignoreBoundaryFaces to ignore cyclics (cyclicACMI)
Remove pyramid edge points, face-diagonals.
bitSet ignoreBoundaryFaces_
Optional boundary faces to ignore.
void size(const label n)
Older name for setAddressableSize.
unsigned int count(const bool on=true) const
Count number of bits set.
isoSurfaceTopo(const polyMesh &mesh, const scalarField &cellValues, const scalarField &pointValues, const scalar iso, const isoSurfaceParams ¶ms=isoSurfaceParams(), const bitSet &ignoreCells=bitSet())
Construct from cell and point values.
#define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2)
label blockCells(UList< cutType > &cuts, const bitSet &ignoreCells) const
Mark ignoreCells as BLOCKED.
Marching tet iso surface algorithm with optional filtering to keep only points originating from mesh ...
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
void resize(const label len)
Adjust allocated size of list.
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
label nPoints() const noexcept
Number of mesh points.
const scalarField & cVals_
Cell values.
void append(const T &val)
Append an element at the end of the list.
virtual const labelList & faceNeighbour() const
Return face neighbour.
List< edge > edgeList
List of edge.
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
Preferences for controlling iso-surface algorithms.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
constexpr char nl
The newline '\n' character (0x0a)
T & first()
Access first element of the list, position [0].
Ostream & endl(Ostream &os)
Add newline and flush stream.
const cellList & cells() const
SubList< label > subList
Declare type of subList.
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
label nFaces() const noexcept
Number of mesh faces.
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
UList< bool > boolUList
A UList of bools.
const labelListList & edgeLoops() const
Return list of closed loops of boundary vertices.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
tmp< Field< Type > > interpolateTemplate(const Field< Type > &cellData, const Field< Type > &pointData) const
Interpolates cellData and pointData fields.
UList< label > labelUList
A UList of labels.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Low-level components common to various iso-surface algorithms.
const scalar iso_
Iso value.
MinMax< label > minMax(const labelHashSet &set)
Find the min/max values of labelHashSet.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
void inplaceSubsetMesh(const bitSet &include)
Subset the surface using the selected faces.
List< face > faceList
List of faces.
pointField vertices(const blockVertexList &bvl)
#define defineIsoSurfaceInterpolateMethods(ThisClass)
const boundBox & getClipBounds() const noexcept
Get optional clipping bounding box.
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
const polyMesh & mesh_
Reference to mesh.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
fileName globalPath() const
Return global path for the case = rootPath/globalCaseName. Same as TimePaths::globalPath() ...
A class for handling words, derived from Foam::string.
const Time & time() const noexcept
Return time registry.
static void appendTriLabels(DynamicList< label > &verts, const label a, const label b, const label c, const bool flip)
label calcCellCuts(List< cutType > &cuts) const
Populate a list of candidate cell cuts using getCellCutType()
static labelList adjustTetBasePtIs(const polyMesh &mesh, const bool report=false)
Return an adjusted list of tet base points.
scalar mag() const
The magnitude/length of the bounding box diagonal.
#define SNAP_END_VALUE(pos, val)
labelList meshCells_
For every face, the original cell in mesh.
virtual const labelList & faceOwner() const
Return face owner.
const polyMesh & mesh() const noexcept
The mesh for which the iso-surface is associated.
const labelListList & edgeFaces() const
Return edge-face addressing.
label nInternalFaces() const noexcept
Number of internal faces.
label timeIndex() const noexcept
Return the current time index.
const Field< point_type > & points() const noexcept
Return reference to global points.
virtual const faceList & faces() const
Return raw faces.
A location outside the volume.
const vectorField & cellCentres() const
bool test(const label pos) const
Test for True value at specified position, never auto-vivify entries.
static word printf(const char *fmt, const PrimitiveType &val)
Use a printf-style formatter for a primitive.
const scalarField & pVals_
Point values.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
T & emplace_back(Args &&... args)
Construct an element at the end of the list, return reference to the new list element.
static int getTetCutIndex(scalar p0, scalar p1, scalar p2, scalar p3, const scalar val, const bool doSnap) noexcept
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
static const Enum< filterType > filterNames
Names for the filtering types.
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field values.
label size() const
The surface size is the number of faces.
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label nCells() const noexcept
Number of mesh cells.
static constexpr Foam::label BLOCKED
const dimensionedScalar c
Speed of light in a vacuum.
List< surfZone > surfZoneList
List of surfZone.
void clearStorage()
Clear the list and delete storage.
static const point & getMeshPointRef(const polyMesh &mesh, const label pointi)
messageStream Info
Information stream (stdout output on master, null elsewhere)
const boundBox & bounds() const noexcept
Return mesh bounding box.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
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))
bool snap() const noexcept
Get point snapping flag.
Convenience macros for instantiating iso-surface interpolate methods.
filterType filter() const noexcept
Get current filter type.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Do not request registration (bool: false)
const volScalarField & p0
const dimensionedScalar mp
Proton mass.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
static constexpr const zero Zero
Global zero (0)