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_.resize(4*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
804 compactFaces.
append(face(loop.size()));
805 face&
f = compactFaces.
last();
810 const label pointi =
mp[loop[i]];
811 if (filterDiag && pointFromDiag[pointi])
813 const label prevPointi =
mp[loop[loop.fcIndex(i)]];
816 pointFromDiag[prevPointi]
817 && (pointToFace[pointi] != pointToFace[prevPointi])
842 const label pointi =
mp[loop[i]];
846 compactCellIDs.
append(cellID);
852 void Foam::isoSurfaceTopo::removeInsidePoints
855 const bool filterDiag,
863 DynamicList<label>& compactCellIDs
868 compactCellIDs.clear();
869 compactCellIDs.reserve(
s.size()/4);
871 DynamicList<face> compactFaces(
s.size()/4);
873 for (label celli = 0; celli < start.size()-1; ++celli)
877 const label nTris = start[celli+1]-start[celli];
883 SubList<face>(
s, nTris, start[celli]),
901 s.swapFaces(compactFaces);
927 Pout<<
"isoSurfaceTopo:" <<
nl 930 <<
" isoValue : " << iso <<
nl 934 <<
" ignoreCells : " << ignoreCells.
count()
941 label nBlockedCells = 0;
944 nBlockedCells +=
blockCells(cellCutType_, ignoreCells);
962 Pout<<
"isoSurfaceTopo : candidate cells cut " 964 <<
" blocked " << nBlockedCells
970 const auto& fvmesh = refCast<const fvMesh>(
mesh);
976 "isoSurfaceTopo.cutType",
977 fvmesh.time().timeName(),
989 forAll(cellCutType_, celli)
991 debugFld[celli] = cellCutType_[celli];
994 Info<<
"Writing cut types: " << debugField.objectRelPath() <<
nl;
1003 constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
1007 forAll(cellCutType_, celli)
1009 if ((cellCutType_[celli] & realCut) != 0)
1011 debugCutCells[nout] = celli;
1013 if (nout >= nCutCells)
break;
1018 vtk::vtuCells vtuCells;
1019 vtuCells.reset(
mesh_, debugCutCells);
1021 vtk::internalMeshWriter
writer 1028 / (
"isoSurfaceTopo." + timeDesc +
"-cutCells")
1043 Info<<
"isoSurfaceTopo : (debug) wrote " 1051 tetCutAddressing tetCutAddr
1059 for (label celli = 0; celli <
mesh_.
nCells(); ++celli)
1061 startTri[celli] = tetCutAddr.nFaces();
1062 if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
1068 bool(cellCutType_[celli] & cutType::TETCUT),
1075 startTri.last() = tetCutAddr.nFaces();
1078 tetBasePtIs.clear();
1079 tetCutAddr.clearHashes();
1083 faceList allTriFaces(startTri.last());
1085 auto& verts = tetCutAddr.cutPoints();
1088 for (face&
f : allTriFaces)
1091 f[0] = verts[verti++];
1092 f[1] = verts[verti++];
1093 f[2] = verts[verti++];
1101 for (label celli = 0; celli < startTri.size()-1; ++celli)
1107 (startTri[celli+1] - startTri[celli]),
1113 pointToVerts_.
transfer(tetCutAddr.pointToVerts());
1128 std::move(allTriPoints),
1129 std::move(allTriFaces),
1135 Pout<<
"isoSurfaceTopo : generated " 1143 const Mesh&
s = *
this;
1145 vtk::surfaceWriter
writer 1152 / (
"isoSurfaceTopo." + timeDesc +
"-triangles")
1173 const edge& verts = pointToVerts_[i];
1174 if (verts.first() == verts.second())
1179 if (tetCutAddr.pointFromDiag().test(i))
1182 pointStatus[i] = -1;
1186 writer.write(
"point-status", pointStatus);
1189 Info<<
"isoSurfaceTopo : (debug) wrote " 1212 s.compactPoints(pointMap);
1213 pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
1224 DynamicList<label> compactCellIDs;
1234 tetCutAddr.pointFromDiag(),
1235 tetCutAddr.pointToFace(),
1241 s.compactPoints(pointMap);
1243 pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
1248 Pout<<
"isoSurfaceTopo :" 1249 " after removing cell centre and face-diag triangles: " 1257 tetCutAddr.clearDiagonal();
1263 bitSet isProtectedPoint;
1267 || tetCutAddr.debugCutTetsOn()
1308 cellCutType_.clear();
1312 if (tetCutAddr.debugCutTetsOn())
1316 const auto& debugCuts = tetCutAddr.debugCutTets();
1319 vtk::vtuCells vtuCells;
1320 vtuCells.resetShapes(debugCuts);
1326 vtk::internalMeshWriter
writer 1333 / (
"isoSurfaceTopo." + timeDesc +
"-cutTets")
1345 Field<scalar> cutTetQuality(debugCuts.size());
1346 forAll(cutTetQuality, teti)
1356 writer.writeCellData(
"tetQuality", cutTetQuality);
1366 for (
const edge& verts : pointToVerts_)
1368 if (verts.first() == verts.second())
1371 pointStatus[verts.first()] = 1;
1375 writer.writePointData(
"point-status", pointStatus);
1378 Info<<
"isoSurfaceTopo : (debug) wrote " 1389 const Mesh&
s = *
this;
1394 openEdgeIds.resize(2*
s.size());
1399 if (eFaces.size() == 1)
1403 const edge&
e = surfEdges[edgei];
1404 const edge& verts0 = pointToVerts_[
mp[
e.first()]];
1405 const edge& verts1 = pointToVerts_[
mp[
e.second()]];
1409 isProtectedPoint.test(verts0.first())
1410 && isProtectedPoint.test(verts0.second())
1411 && isProtectedPoint.test(verts1.first())
1412 && isProtectedPoint.test(verts1.second())
1420 openEdgeIds.insert(edgei);
1425 const label nOpenEdges
1435 openEdgeIds.sortedToc()
1445 / (
"isoSurfaceTopo." + timeDesc +
"-openEdges")
1455 Info<<
"isoSurfaceTopo : (debug) wrote " 1456 << nOpenEdges <<
" open edges: " 1461 Info<<
"isoSurfaceTopo : no open edges" <<
nl;
1468 const Mesh& s = *
this;
1470 vtk::surfaceWriter
writer 1477 / (
"isoSurfaceTopo." + timeDesc +
"-surface")
1498 const edge& verts = pointToVerts_[i];
1499 if (verts.first() == verts.second())
1506 writer.write(
"point-status", pointStatus);
1509 Info<<
"isoSurfaceTopo : (debug) wrote " 1515 tetCutAddr.clearDebug();
1532 bitSet faceSelection;
1539 UIndirectList<face>(surf, faceAddr),
1543 faceSelection.clear();
1544 faceSelection.resize(erosion.size());
1547 const edgeList& surfEdges = erosion.edges();
1550 label nEdgeRemove = 0;
1555 if (eFaces.size() == 1)
1559 const edge&
e = surfEdges[edgei];
1560 const edge& verts0 = pointToVerts_[
mp[
e.first()]];
1561 const edge& verts1 = pointToVerts_[
mp[
e.second()]];
1565 isProtectedPoint.test(verts0.first())
1566 && isProtectedPoint.test(verts0.second())
1567 && isProtectedPoint.test(verts1.first())
1568 && isProtectedPoint.test(verts1.second())
1576 faceSelection.set(eFaces[0]);
1584 Pout<<
"isoSurfaceTopo :" 1585 <<
" removing " << faceSelection.count()
1586 <<
" / " << faceSelection.size()
1587 <<
" faces on " << nEdgeRemove <<
" open edges" <<
endl;
1603 if (surf.size() != faceAddr.size())
1605 faceSelection.clear();
1606 faceSelection.resize(surf.size());
1607 faceSelection.set(faceAddr);
1623 Mesh::subsetMesh(include, pointMap,
faceMap)
1625 Mesh::transfer(filtered);
1627 meshCells_ = UIndirectList<label>(meshCells_,
faceMap)();
1629 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)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
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 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 & last()
Access last element of the list, position [size()-1].
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())
fileName globalPath() const
Return global path for the case.
static constexpr const zero Zero
Global zero (0)