44 const scalar surfaceFeatures::parallelTolerance =
sin(
degToRad(1.0));
76 mag(eHit.rawPoint() - start)
77 <
mag(eHit.rawPoint() - end)
99 for (
label i = 0; i < externalStart_; i++)
101 edgeStat[featureEdges_[i]] = REGION;
105 for (
label i = externalStart_; i < internalStart_; i++)
107 edgeStat[featureEdges_[i]] = EXTERNAL;
111 for (
label i = internalStart_; i < featureEdges_.size(); i++)
113 edgeStat[featureEdges_[i]] = INTERNAL;
124 const scalar includedAngle
135 if (edgeStat[edgeI] == REGION)
139 else if (edgeStat[edgeI] == EXTERNAL)
143 else if (edgeStat[edgeI] == INTERNAL)
149 externalStart_ = nRegion;
150 internalStart_ = externalStart_ + nExternal;
155 featureEdges_.setSize(internalStart_ + nInternal);
158 label externalI = externalStart_;
159 label internalI = internalStart_;
163 if (edgeStat[edgeI] == REGION)
165 featureEdges_[regionI++] = edgeI;
167 else if (edgeStat[edgeI] == EXTERNAL)
169 featureEdges_[externalI++] = edgeI;
171 else if (edgeStat[edgeI] == INTERNAL)
173 featureEdges_[internalI++] = edgeI;
179 calcFeatPoints(edgeStat, minCos);
184 void Foam::surfaceFeatures::calcFeatPoints
193 const edgeList& edges = surf_.edges();
194 const pointField& localPoints = surf_.localPoints();
196 forAll(pointEdges, pointI)
198 const labelList& pEdges = pointEdges[pointI];
200 label nFeatEdges = 0;
204 if (edgeStat[pEdges[i]] != NONE)
212 featurePoints.
append(pointI);
214 else if (nFeatEdges == 2)
221 const label edgeI = pEdges[i];
223 if (edgeStat[edgeI] != NONE)
225 edgeVecs.
append(edges[edgeI].vec(localPoints));
230 if (
mag(edgeVecs[0] & edgeVecs[1]) < minCos)
232 featurePoints.
append(pointI);
237 featurePoints_.transfer(featurePoints);
241 void Foam::surfaceFeatures::classifyFeatureAngles
246 const bool geometricTestOnly
249 const vectorField& faceNormals = surf_.faceNormals();
253 bool selectAll = (
mag(minCos-1.0) < SMALL);
257 const labelList& eFaces = edgeFaces[edgeI];
259 if (eFaces.
size() != 2)
262 edgeStat[edgeI] = REGION;
266 label face0 = eFaces[0];
267 label face1 = eFaces[1];
272 && surf_[face0].region() != surf_[face1].region()
275 edgeStat[edgeI] = REGION;
280 || ((faceNormals[face0] & faceNormals[face1]) < minCos)
286 surf_[face1].
centre(points)
287 - surf_[face0].
centre(points);
289 if ((f0Tof1 & faceNormals[face0]) >= 0.0)
291 edgeStat[edgeI] = INTERNAL;
295 edgeStat[edgeI] = EXTERNAL;
308 const label unsetVal,
309 const label prevEdgeI,
313 const labelList& pEdges = surf_.pointEdges()[vertI];
315 label nextEdgeI = -1;
319 label edgeI = pEdges[i];
324 && edgeStat[edgeI] != NONE
325 && featVisited[edgeI] == unsetVal
352 Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
356 const label startEdgeI,
357 const label startPointI,
358 const label currentFeatI,
362 label edgeI = startEdgeI;
364 label vertI = startPointI;
366 scalar visitedLength = 0.0;
370 if (
findIndex(featurePoints_, startPointI) >= 0)
374 return labelScalar(nVisited, visitedLength);
394 unsetVal = currentFeatI;
400 edgeI = nextFeatEdge(edgeStat, featVisited, unsetVal, edgeI, vertI);
402 if (edgeI == -1 || edgeI == startEdgeI)
411 featVisited[edgeI] = currentFeatI;
415 featVisited[edgeI] = -2;
420 const edge&
e = surf_.edges()[edgeI];
426 visitedLength += e.
mag(surf_.localPoints());
430 if (nVisited > surf_.nEdges())
432 Warning<<
"walkSegment : reached iteration limit in walking " 433 <<
"feature edges on surface from edge:" << startEdgeI
434 <<
" vertex:" << startPointI <<
nl 435 <<
"Returning with large length" <<
endl;
437 return labelScalar(nVisited, GREAT);
442 return labelScalar(nVisited, visitedLength);
464 const label externalStart,
465 const label internalStart
469 featurePoints_(featurePoints),
470 featureEdges_(featureEdges),
471 externalStart_(externalStart),
472 internalStart_(externalStart)
480 const scalar includedAngle,
482 const label minElems,
483 const bool geometricTestOnly
494 if (minLen > 0 || minElems > 0)
509 featurePoints_(featInfoDict.
lookup(
"featurePoints")),
510 featureEdges_(featInfoDict.
lookup(
"featureEdges")),
545 const scalar mergeTol,
546 const bool geometricTestOnly
560 scalar mergeTolSqr =
sqr(mergeTol);
578 const label sEdge = edgeLabel[sEdgeI];
585 dynFeatEdges.
insert(surfEdges[sEdge], count++);
586 dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
592 classifyFeatureAngles
608 if (iter != dynFeatEdges.
end())
610 allEdgeStat[eI] = edgeStat[iter()];
615 dynFeatEdges.clear();
625 featurePoints_(sf.featurePoints()),
626 featureEdges_(sf.featureEdges()),
627 externalStart_(sf.externalStart()),
628 internalStart_(sf.internalStart())
636 const bool regionEdges,
637 const bool externalEdges,
638 const bool internalEdges
647 for (
label i = 0; i < externalStart_; i++)
649 selectedEdges.
append(featureEdges_[i]);
657 for (
label i = externalStart_; i < internalStart_; i++)
659 selectedEdges.
append(featureEdges_[i]);
667 for (
label i = internalStart_; i < featureEdges_.
size(); i++)
669 selectedEdges.
append(featureEdges_[i]);
673 return selectedEdges.
shrink();
679 const scalar includedAngle,
680 const bool geometricTestOnly
688 classifyFeatureAngles
705 const label minElems,
706 const scalar includedAngle
723 label startEdgeI = 0;
728 for (; startEdgeI < edgeStat.
size(); startEdgeI++)
732 edgeStat[startEdgeI] !=
NONE 733 && featLines[startEdgeI] == -1
741 if (startEdgeI == edgeStat.
size())
748 featLines[startEdgeI] = featI;
750 const edge& startEdge = surf_.
edges()[startEdgeI];
753 labelScalar leftPath =
764 labelScalar rightPath =
783 || (leftPath.n_ + rightPath.n_ + 1 < minElems)
789 featLines[startEdgeI] = -2;
821 label edgeI = featureEdges_[i];
823 if (featLines[edgeI] == -2)
825 edgeStat[edgeI] =
NONE;
840 featInfoDict.
add(
"externalStart", externalStart_);
841 featInfoDict.
add(
"internalStart", internalStart_);
842 featInfoDict.
add(
"featureEdges", featureEdges_);
843 featInfoDict.
add(
"featurePoints", featurePoints_);
845 featInfoDict.
write(writeFile);
859 OFstream regionStr(prefix +
"_regionEdges.obj");
860 Pout<<
"Writing region edges to " << regionStr.
name() <<
endl;
863 for (
label i = 0; i < externalStart_; i++)
865 const edge&
e = surf_.
edges()[featureEdges_[i]];
869 regionStr <<
"l " << verti-1 <<
' ' << verti <<
endl;
873 OFstream externalStr(prefix +
"_externalEdges.obj");
874 Pout<<
"Writing external edges to " << externalStr.
name() <<
endl;
877 for (
label i = externalStart_; i < internalStart_; i++)
879 const edge&
e = surf_.
edges()[featureEdges_[i]];
883 externalStr <<
"l " << verti-1 <<
' ' << verti <<
endl;
886 OFstream internalStr(prefix +
"_internalEdges.obj");
887 Pout<<
"Writing internal edges to " << internalStr.
name() <<
endl;
890 for (
label i = internalStart_; i < featureEdges_.
size(); i++)
892 const edge&
e = surf_.
edges()[featureEdges_[i]];
896 internalStr <<
"l " << verti-1 <<
' ' << verti <<
endl;
899 OFstream pointStr(prefix +
"_points.obj");
900 Pout<<
"Writing feature points to " << pointStr.
name() <<
endl;
904 label pointI = featurePoints_[i];
940 label surfPointI = pointLabels[i];
942 const point& surfPt = surfPoints[surfPointI];
953 <<
"Problem for point " 954 << surfPointI <<
" in tree " << ppTree.
bb()
960 if (
magSqr(samples[sampleI] - surfPt) < maxDistSqr[sampleI])
962 nearest.insert(sampleI, surfPointI);
974 <<
"Dumping nearest surface feature points to nearestSamples.obj" 976 <<
"View this Lightwave-OBJ file with e.g. javaview" <<
endl 979 OFstream objStream(
"nearestSamples.obj");
986 objStream<<
"l " << vertI-1 <<
' ' << vertI <<
endl;
1002 const scalar minSampleDist
1008 scalar maxSearchSqr =
max(maxDistSqr);
1027 label surfEdgeI = selectedEdges[i];
1029 const edge&
e = surfEdges[surfEdgeI];
1031 if (debug && (i % 1000) == 0)
1033 Pout<<
"looking at surface feature edge " << surfEdgeI
1034 <<
" verts:" << e <<
" points:" << surfPoints[e[0]]
1035 <<
' ' << surfPoints[e[1]] <<
endl;
1040 scalar eMag =
mag(eVec);
1055 point edgePoint(surfPoints[e.
start()] + s*eVec);
1069 label sampleI = info.index();
1071 if (
magSqr(info.hitPoint() - edgePoint) < maxDistSqr[sampleI])
1073 nearest.insert(sampleI, surfEdgeI);
1083 s +=
max(minSampleDist*eMag, sampleDist[sampleI]);
1085 if (s >= (1-minSampleDist)*eMag)
1100 Pout<<
"Dumping nearest surface edges to nearestEdges.obj\n" 1101 <<
"View this Lightwave-OBJ file with e.g. javaview\n" <<
endl;
1103 OFstream objStream(
"nearestEdges.obj");
1108 const label sampleI = iter.key();
1112 const edge&
e = surfEdges[iter()];
1119 objStream<<
"l " << vertI-1 <<
' ' << vertI <<
endl;
1140 const scalar minSampleDist
1162 scalar maxSearchSqr =
max(maxDistSqr);
1173 label surfEdgeI = selectedEdges[i];
1175 const edge&
e = surfEdges[surfEdgeI];
1177 if (debug && (i % 1000) == 0)
1179 Pout<<
"looking at surface feature edge " << surfEdgeI
1180 <<
" verts:" << e <<
" points:" << surfPoints[e[0]]
1181 <<
' ' << surfPoints[e[1]] <<
endl;
1186 scalar eMag =
mag(eVec);
1201 point edgePoint(surfPoints[e.
start()] + s*eVec);
1215 label index = info.index();
1217 label sampleEdgeI = ppTree.
shapes().edgeLabels()[index];
1219 const edge& e = sampleEdges[sampleEdgeI];
1221 if (
magSqr(info.hitPoint() - edgePoint) < maxDistSqr[e.
start()])
1240 if (s >= (1-minSampleDist)*eMag)
1254 Pout<<
"Dumping nearest surface feature edges to nearestEdges.obj\n" 1255 <<
"View this Lightwave-OBJ file with e.g. javaview\n" <<
endl;
1257 OFstream objStream(
"nearestEdges.obj");
1262 const label sampleEdgeI = iter.key();
1264 const edge& sampleEdge = sampleEdges[sampleEdgeI];
1273 objStream<<
"l " << vertI-1 <<
' ' << vertI <<
endl;
1287 scalar searchSpanSqr,
1319 const point& sample = samples[i];
1333 edgeLabel[i] = selectedEdges[info.
index()];
1340 localPoints[e.
start()],
1341 localPoints[e.
end()],
1346 edgeEndPoint[i] = pHit.
index();
1361 const vector& searchSpan,
1369 pointOnFeature.
setSize(selectedSampleEdges.
size());
1388 forAll(selectedSampleEdges, i)
1390 const edge&
e = sampleEdges[selectedSampleEdges[i]];
1396 treeBoundBox tightest(eMid - searchSpan, eMid + searchSpan);
1411 edgeLabel[i] = selectedEdges[info.
index()];
1413 pointOnFeature[i] = info.
hitPoint();
1423 scalar searchSpanSqr,
1430 searchDomain.inflate(0.1);
1452 const edge& sample = surfEdges[edgeI];
1454 const point& startPoint = surfLocalPoints[sample.
start()];
1465 const vector surfEdgeDir = midPoint - startPoint;
1467 const edge& featEdge = edges[infoMid.
index()];
1468 const vector featEdgeDir = featEdge.
vec(points);
1471 if (
mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
1473 edgeLabel[edgeI] = edgeI;
1489 "Foam::surfaceFeatures::operator=(const Foam::surfaceFeatures&)" 1490 ) <<
"Attempted assignment to self" 1498 "Foam::surfaceFeatures::operator=(const Foam::surfaceFeatures&)" 1499 ) <<
"Operating on different surfaces"
PointHit< Point > nearestDist(const Point &p) const
Return nearest distance to line from a given point.
void setCapacity(const label)
Alter the size of the underlying storage.
void write(const fileName &fName) const
Write as dictionary to file.
const labelList & featureEdges() const
Return feature edge list.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
vector point
Point is a vector.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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 ))
const labelList & featurePoints() const
Return feature point list.
dimensioned< scalar > mag(const dimensioned< Type > &)
void nearestFeatEdge(const edgeList &edges, const pointField &points, scalar searchSpanSqr, labelList &edgeLabel) const
Find nearest feature edge to each surface edge. Uses the.
labelList selectFeatureEdges(const bool regionEdges, const bool externalEdges, const bool internalEdges) const
Helper function: select a subset of featureEdges_.
const labelListList & edgeFaces() const
Return edge-face addressing.
label externalStart() const
Start of external edges.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
bool insert(const edge &, const T &newElmt)
Insert a new hashedEntry.
const Type & shapes() const
Reference to shape.
const Point & rawPoint() const
Return point with no checking.
T & last()
Return the last element of the list.
A HashTable to objects of type <T> with a label key.
const treeBoundBox & bb() const
Top bounding box.
Triangulated surface description with patch information.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Standard boundBox + extra functionality for use in octree.
void size(const label)
Override size to be inconsistent with allocated storage.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
surfaceFeatures(const triSurface &)
Construct from surface.
Map< pointIndexHit > nearestEdges(const labelList &selectedEdges, const edgeList &sampleEdges, const labelList &selectedSampleEdges, const pointField &samplePoints, const scalarField &sampleDist, const scalarField &maxDistSqr, const scalar minSampleDist=0.1) const
Like nearestSamples but now gets nearest point on.
PointIndexHit< point > pointIndexHit
bool hit() const
Is there a hit.
A list of keyword definitions, which are a keyword followed by any number of values (e...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Holds feature edges/points of surface.
label index() const
Return index.
volScalarField sf(fieldObject, mesh)
PointHit< point > pointHit
label readLabel(Istream &is)
const Point & hitPoint() const
Return hit point.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
scalar mag(const pointField &) const
Return scalar magnitude.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
const Field< PointType > & localPoints() const
Return pointField of points in patch.
Map< label > nearestSamples(const labelList &selectedPoints, const pointField &samples, const scalarField &maxDistSqr) const
Find nearest sample for selected surface points.
void clear()
Clear the list, i.e. set size to zero.
const double e
Elementary charge.
void setSize(const label)
Reset size of List.
Ostream & endl(Ostream &os)
Add newline and flush stream.
linePointRef line(const pointField &) const
Return edge line.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Mid-point interpolation (weighting factors = 0.5) scheme class.
const fileName & name() const
Return the name of the stream.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
labelList trimFeatures(const scalar minLen, const label minElems, const scalar includedAngle)
Delete small sets of edges. Edges are stringed up and any.
dimensionedScalar cos(const dimensionedScalar &ds)
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
label otherVertex(const label a) const
Given one vertex, return the other.
void nearestSurfEdge(const labelList &selectedEdges, const pointField &samples, scalar searchSpanSqr, labelList &edgeLabel, labelList &edgeEndPoint, pointField &edgePoint) const
Find nearest surface edge (out of selectedEdges) for.
Unit conversion functions.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
void findFeatures(const scalar includedAngle, const bool geometricTestOnly)
Find feature edges using provided included angle.
bool add(entry *, bool mergeEntry=false)
Add a new entry.
vector vec(const pointField &) const
Return the vector (end - start)
label nRegionEdges() const
Return number of region edges.
errorManip< error > abort(error &err)
void writeDict(Ostream &) const
Write as dictionary.
const triSurface & surface() const
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt > > &) const
Return *this (used for point which is a typedef to Vector<scalar>.
List< edgeStatus > toStatus() const
From member feature edges to status per edge.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Non-pointer based hierarchical recursive searching.
Holds data for octree to work on an edges subset.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
label end() const
Return end vertex label.
List< label > labelList
A List of labels.
A class for handling file names.
label nExternalEdges() const
Return number of external edges.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Point centre() const
Return centre (centroid)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
void append(const T &)
Append an element at the end of the list.
label start() const
Return start vertex label.
line< point, const point & > linePointRef
Line using referred points.
void writeObj(const fileName &prefix) const
Write to separate OBJ files (region, external, internal edges,.
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
label nInternalEdges() const
Return number of internal edges.
void write(Ostream &, const bool subDict=true) const
Write dictionary, normally with sub-dictionary formatting.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void setFromStatus(const List< edgeStatus > &edgeStat, const scalar includedAngle)
Set from status per edge.
label nEdges() const
Return number of edges in patch.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Holds (reference to) pointField. Encapsulation of data needed for octree searches. Used for searching for nearest point. No bounding boxes around points. Only overlaps and calcNearest are implemented, rest makes little sense.
void operator=(const surfaceFeatures &)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
label internalStart() const
Start of internal edges.
point centre(const pointField &) const
Return centre (centroid)
dimensionedScalar sin(const dimensionedScalar &ds)