45 {
"none",
"partial",
"full"};
53 bool Foam::isoSurface::isTriCut
59 bool aLower = (pointValues[tri[0]] < iso_);
60 bool bLower = (pointValues[tri[1]] < iso_);
61 bool cLower = (pointValues[tri[2]] < iso_);
63 return !(aLower == bLower && aLower == cLower);
67 Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
73 const cell& cFaces = mesh_.cells()[celli];
79 const face&
f = mesh_.faces()[cFaces[cFacei]];
81 for (
label fp = 1; fp < f.
size() - 1; fp++)
85 if (isTriCut(tri, pVals_))
87 return cellCutType::cut;
91 return cellCutType::notCut;
95 bool cellLower = (cVals_[celli] < iso_);
102 label facei = cFaces[cFacei];
103 const face&
f = mesh_.faces()[facei];
108 if ((pVals_[f[fp]] < iso_) != cellLower)
120 const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
127 if (isTriCut(
triFace(f[fp0], f[fp], f[nextFp]), pVals_))
148 const labelList& cPoints = mesh_.cellPoints(celli);
154 if ((pVals_[cPoints[i]] < iso_) != cellLower)
160 if (nPyrCuts == cPoints.
size())
162 return cellCutType::sphere;
166 return cellCutType::cut;
171 return cellCutType::notCut;
189 cellCutTypes[celli] = calcCutType(tet.
isA(mesh_, celli), celli);
191 if (cellCutTypes[celli] == cellCutType::cut)
199 Pout<<
"isoSurface : detected " << nCutCells
200 <<
" candidate cut cells." <<
endl;
206 Foam::scalar Foam::isoSurface::minTetQ
209 const label faceBasePtI
215 mesh_.cellCentres()[mesh_.faceOwner()[facei]],
221 if (mesh_.isInternalFace(facei))
229 mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
240 void Foam::isoSurface::fixTetBasePtIs()
244 const faceList& faces = mesh_.faces();
245 const labelList& faceOwner = mesh_.faceOwner();
246 const labelList& faceNeighbour = mesh_.faceNeighbour();
250 tetBasePtIs_ = mesh_.tetBasePtIs();
255 forAll(tetBasePtIs_, facei)
257 if (tetBasePtIs_[facei] == -1)
259 problemCells[faceOwner[facei]] =
true;
260 if (mesh_.isInternalFace(facei))
262 problemCells[faceNeighbour[facei]] =
true;
273 if (problemCells[celli])
275 const cell& cFaces = cells[celli];
280 const label facei = cFaces[i];
281 const face&
f = faces[facei];
284 const label pointi = f[fp];
287 if (pointFnd == pointCount.
end())
289 pointCount.
insert(pointi, 1);
303 <<
" at:" << mesh_.points()[iter.key()]
308 problemPoints[iter.key()] =
true;
319 forAll(tetBasePtIs_, facei)
323 problemCells[faceOwner[facei]]
324 || (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
327 const face&
f = faces[facei];
331 const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
333 const bool prevPointIsProblem = problemPoints[f[f.
rcIndex(fp0)]];
334 const bool nextPointIsProblem = problemPoints[f[f.fcIndex(fp0)]];
336 if (!prevPointIsProblem && !nextPointIsProblem)
344 scalar maxQ = -
GREAT;
348 const bool prevPointIsProblem = problemPoints[f[f.rcIndex(fp)]];
349 const bool nextPointIsProblem = problemPoints[f[f.fcIndex(fp)]];
351 if (!prevPointIsProblem && !nextPointIsProblem)
353 const scalar q = minTetQ(facei, fp);
365 tetBasePtIs_[facei] = maxFp;
380 Pout<<
"isoSurface : adapted starting point of triangulation on " 381 << nAdapted <<
" faces." <<
endl;
389 const bool edgeIsDiag,
390 const edge& vertices,
399 if (edgeFnd != vertsToPoint.
end())
408 pointToVerts.
append(vertices);
409 pointToFace.
append(facei);
410 pointFromDiag.
append(edgeIsDiag);
411 vertsToPoint.
insert(vertices, pointi);
417 void Foam::isoSurface::generateTriPoints
467 edge(pIndex[0], pIndex[1]),
468 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
477 edge(pIndex[0], pIndex[2]),
478 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
487 edge(pIndex[0], pIndex[3]),
488 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
492 if (triIndex == 0x0E)
496 Swap(verts[sz-2], verts[sz-1]);
510 edge(pIndex[1], pIndex[0]),
511 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
520 edge(pIndex[1], pIndex[3]),
521 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
530 edge(pIndex[1], pIndex[2]),
531 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
535 if (triIndex == 0x0D)
539 Swap(verts[sz-2], verts[sz-1]);
553 edge(pIndex[0], pIndex[2]),
554 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
563 edge(pIndex[1], pIndex[3]),
564 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
574 edge(pIndex[0], pIndex[3]),
575 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
587 edge(pIndex[1], pIndex[2]),
588 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
593 if (triIndex == 0x0C)
597 Swap(verts[sz-5], verts[sz-4]);
598 Swap(verts[sz-2], verts[sz-1]);
612 edge(pIndex[2], pIndex[0]),
613 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
622 edge(pIndex[2], pIndex[1]),
623 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
632 edge(pIndex[2], pIndex[3]),
633 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
637 if (triIndex == 0x0B)
641 Swap(verts[sz-2], verts[sz-1]);
655 edge(pIndex[0], pIndex[1]),
656 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
665 edge(pIndex[2], pIndex[3]),
666 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
678 edge(pIndex[0], pIndex[3]),
679 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
689 edge(pIndex[1], pIndex[2]),
690 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
695 if (triIndex == 0x0A)
699 Swap(verts[sz-5], verts[sz-4]);
700 Swap(verts[sz-2], verts[sz-1]);
714 edge(pIndex[0], pIndex[1]),
715 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
724 edge(pIndex[2], pIndex[3]),
725 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
736 edge(pIndex[1], pIndex[3]),
737 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
749 edge(pIndex[0], pIndex[2]),
750 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
754 if (triIndex == 0x09)
758 Swap(verts[sz-5], verts[sz-4]);
759 Swap(verts[sz-2], verts[sz-1]);
773 edge(pIndex[3], pIndex[0]),
774 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
783 edge(pIndex[3], pIndex[2]),
784 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
793 edge(pIndex[3], pIndex[1]),
794 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
797 if (triIndex == 0x07)
801 Swap(verts[sz-2], verts[sz-1]);
809 void Foam::isoSurface::generateTriPoints
823 const faceList& faces = mesh_.faces();
824 const labelList& faceOwner = mesh_.faceOwner();
826 const cell& cFaces = mesh_.cells()[celli];
833 label facei = cFaces[0];
834 const face& f0 = faces[facei];
837 const face&
f1 = faces[cFaces[1]];
838 label oppositeI = -1;
854 if (faceOwner[facei] == celli)
901 label nTris = (verts.
size()-startTrii)/3;
902 for (
label i = 0; i < nTris; i++)
909 const pointField& cellCentres = mesh_.cellCentres();
913 label facei = cFaces[cFacei];
914 const face&
f = faces[facei];
916 const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
929 label p2 = f[nextFp];
930 if (faceOwner[facei] == celli)
933 if (i != 2) edgeIsDiag[1] =
true;
934 if (i != f.
size()-1) edgeIsDiag[0] =
true;
938 if (i != 2) edgeIsDiag[0] =
true;
939 if (i != f.
size()-1) edgeIsDiag[1] =
true;
972 mesh_.nPoints()+celli
986 label nTris = (verts.
size()-startTrii)/3;
987 for (
label i = 0; i < nTris; i++)
996 void Foam::isoSurface::triangulateOutside
998 const bool filterDiag,
1020 const labelList& loop = edgeLoops[loopi];
1022 if (loop.
size() > 2)
1031 label pointi = mp[loop[i]];
1032 if (filterDiag && pointFromDiag[pointi])
1034 label prevPointi = mp[loop[loop.fcIndex(i)]];
1037 pointFromDiag[prevPointi]
1038 && (pointToFace[pointi] != pointToFace[prevPointi])
1063 label pointi = mp[loop[i]];
1067 compactCellIDs.
append(cellID);
1075 const bool filterDiag,
1086 pointCompactMap.
clear();
1087 compactCellIDs.
clear();
1091 for (
label celli = 0; celli < start.
size()-1; celli++)
1095 label nTris = start[celli+1]-start[celli];
1122 pointCompactMap.
clear();
1126 face&
f = compactFaces[i];
1129 label pointi = f[fp];
1130 label compacti = oldToCompact[pointi];
1133 compacti = compactPoints.
size();
1134 oldToCompact[pointi] = compacti;
1135 compactPoints.append(points[pointi]);
1136 pointCompactMap.
append(pointi);
1146 move(compactPoints),
1173 Pout<<
"isoSurface : iso:" << iso_
1174 <<
" filter:" << filterTypeNames_[filter] <<
endl;
1183 const label nCutCells = calcCutTypes(tet, cellCutTypes);
1203 labelList startTri(mesh_.nCells()+1, 0);
1205 for (
label celli = 0; celli < mesh_.nCells(); celli++)
1207 startTri[celli] = faceLabels.
size();
1208 if (cellCutTypes[celli] != cellCutType::notCut)
1213 tet.
isA(mesh_, celli),
1224 for (
label i = startTri[celli]; i < faceLabels.
size(); i++)
1226 cellLabels.
append(celli);
1230 startTri[mesh_.nCells()] = faceLabels.
size();
1233 pointToVerts_.transfer(pointToVerts);
1234 meshCells_.transfer(cellLabels);
1235 pointToFace_.transfer(pointToFace);
1241 mesh_.cellCentres(),
1252 allTris[i].setSize(3);
1253 allTris[i][0] = verts[verti++];
1254 allTris[i][1] = verts[verti++];
1255 allTris[i][2] = verts[verti++];
1270 move(allPoints.
ref()),
1287 Pout<<
"isoSurface : generated " << size() <<
" faces." <<
endl;
1291 if (filter != filterType::none)
1301 (filter == filterType::full ?
true :
false),
1314 meshCells_.transfer(compactCellIDs);
1318 Pout<<
"isoSurface :" 1319 <<
" after removing cell centre and face-diag triangles : " 1324 if (filter == filterType::full)
1345 isBoundaryPoint.
set(mesh.
faces()[facei]);
1360 const labelList& eFaces = edgeFaces[edgei];
1361 if (eFaces.
size() == 1)
1365 const edge&
e = edges()[edgei];
1366 const edge& verts0 = pointToVerts_[mp[e[0]]];
1367 const edge& verts1 = pointToVerts_[mp[e[1]]];
1370 isBoundaryPoint[verts0[0]]
1371 && isBoundaryPoint[verts0[1]]
1372 && isBoundaryPoint[verts1[0]]
1373 && isBoundaryPoint[verts1[1]]
1381 if (removeFace.
set(eFaces[0]))
1392 Pout<<
"isoSurface :" 1393 <<
" removing " << nFaces
1394 <<
" faces since on open edges" <<
endl;
1404 forAll(removeFace, facei)
1406 if (!removeFace[facei])
static const scalar GREAT
A cellMatcher for tet cells.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
A face is a list of labels corresponding to mesh vertices.
A 1D vector of objects of type <T> with a fixed size <Size>.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nInternalFaces() const
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
T & ref() const
Return non-const reference or generate a fatal error.
A surface zone on a MeshedSurface.
void size(const label)
Override size to be inconsistent with allocated storage.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
static scalar minQuality(const polyMesh &mesh, const point &cC, const label fI, const bool isOwner, const label faceBasePtI)
Given a face and cc and starting index for triangulation determine.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const dimensionedScalar & mp
Proton mass.
void transfer(MeshedSurface< Face > &)
Transfer the contents of the argument and annul the argument.
bool insert(const Key &key)
Insert a new entry.
Initialise the NamedEnum HashTable from the static list of names.
const labelListList & edgeLoops() const
Return list of closed loops of boundary vertices.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
A list of faces which address into the list of points.
A List obtained as a section of another List.
void set(const PackedList< 1 > &)
Set specified bits.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A triangular face using a FixedList of labels corresponding to mesh vertices.
static const NamedEnum< filterType, 3 > filterTypeNames_
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
const Field< PointType > & points() const
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
virtual const faceList & faces() const
Return raw faces.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
virtual void reset(pointField &&points, List< Face > &&faces, surfZoneList &&zones)
Reset primitive data (points, faces and zones)
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
List< surfZone > surfZoneList
virtual bool isA(const primitiveMesh &mesh, const label celli)
Exact match. Uses faceSizeMatch.
void setSize(const label)
Reset size of List.
label size() const
The surface size is the number of faces.
A cell is defined as a list of faces with extra functionality.
prefixOSstream Pout(cout, "Pout")
A List with indirect addressing.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const doubleScalar e
Elementary charge.
Mesh consisting of general polyhedral cells.
A class for managing temporary objects.
T & last()
Return the last element of the list.
void clear()
Clear the addressed list, i.e. set the size to zero.
isoSurface(const polyMesh &mesh, const scalarField &cellValues, const scalarField &pointValues, const scalar iso, const filterType filter)
Construct from dictionary.
const List< surfZone > & surfZones() const
Const access to the surface zones.