36 scalar treeDataFace::tolSqr =
sqr(1
e-6);
46 const face& f = mesh_.faces()[facei];
48 treeBoundBox bb(points[f[0]], points[f[0]]);
50 for (
label fp = 1; fp < f.size(); fp++)
52 const point& p = points[f[fp]];
55 bb.max() =
max(bb.max(),
p);
61 void Foam::treeDataFace::update()
65 isTreeFace_.set(faceLabels_[i], 1);
70 bbs_.setSize(faceLabels_.size());
74 bbs_[i] = calcBb(faceLabels_[i]);
90 faceLabels_(faceLabels),
91 isTreeFace_(mesh.
nFaces(), 0),
106 faceLabels_(faceLabels),
107 isTreeFace_(mesh.
nFaces(), 0),
121 faceLabels_(
identity(mesh_.nFaces())),
122 isTreeFace_(mesh.
nFaces(), 0),
141 isTreeFace_(mesh_.nFaces(), 0),
174 cc[i] = mesh_.faceCentres()[faceLabels_[i]];
201 if (info.
index() == -1)
204 <<
"Could not find " << sample <<
" in octree." 214 Pout<<
"getSampleType : sample:" << sample
215 <<
" nearest face:" << facei;
223 const face& f = mesh_.faces()[facei];
224 const vector& area = mesh_.faceAreas()[facei];
225 const point& fc = mesh_.faceCentres()[facei];
240 Pout<<
" -> face hit:" << curPt
241 <<
" comparing to face normal " << area <<
endl;
248 Pout<<
" -> face miss:" << curPt;
256 const scalar typDimSqr =
mag(area) + VSMALL;
260 if ((
magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
272 if (isTreeFace_.get(pFaces[i]) == 1)
274 vector n = mesh_.faceAreas()[pFaces[i]];
275 n /=
mag(n) + VSMALL;
283 Pout<<
" -> face point hit :" << points[f[fp]]
284 <<
" point normal:" << pointNormal
286 <<
magSqr(points[f[fp]] - curPt)/typDimSqr <<
endl;
295 if ((
magSqr(fc - curPt)/typDimSqr) < tolSqr)
302 Pout<<
" -> centre hit:" << fc
303 <<
" distance:" <<
magSqr(fc - curPt)/typDimSqr <<
endl;
315 const labelList& myEdges = mesh_.faceEdges()[facei];
319 const edge&
e = mesh_.edges()[myEdges[myEdgeI]];
326 ).nearestDist(sample);
335 const labelList& eFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
341 if (isTreeFace_.get(eFaces[i]) == 1)
343 vector n = mesh_.faceAreas()[eFaces[i]];
344 n /=
mag(n) + VSMALL;
353 <<
" comparing to edge normal:" << edgeNormal
378 ).nearestDist(sample);
386 vector ePrev = points[f[f.rcIndex(fp)]] - fc;
387 vector eNext = points[f[f.fcIndex(fp)]] - fc;
390 nLeft /=
mag(nLeft) + VSMALL;
392 vector nRight = e ^ eNext;
393 nRight /=
mag(nRight) + VSMALL;
397 Pout<<
" -> internal edge hit point:" << edgeHit.
rawPoint()
398 <<
" comparing to edge normal " 399 << 0.5*(nLeft + nRight)
407 0.5*(nLeft + nRight),
415 Pout<<
"Did not find sample " << sample
416 <<
" anywhere related to nearest face " << facei <<
endl 421 Pout<<
" vertex:" << f[fp] <<
" coord:" << points[f[fp]]
452 if (!cubeBb.
overlaps(calcBb(faceLabels_[index])))
462 label facei = faceLabels_[index];
464 const face& f = mesh_.faces()[facei];
472 const point& fc = mesh_.faceCentres()[facei];
493 void Foam::treeDataFace::findNearestOp::operator()
498 scalar& nearestDistSqr,
507 const label index = indices[i];
514 if (distSqr < nearestDistSqr)
516 nearestDistSqr = distSqr;
524 void Foam::treeDataFace::findNearestOp::operator()
539 bool Foam::treeDataFace::findIntersectOp::operator()
544 point& intersectionPoint
561 const label facei = shape.faceLabels_[index];
563 const vector dir(end - start);
578 intersectionPoint = inter.
hitPoint();
label end() const
Return end vertex label.
const primitiveMesh & mesh() const
A simple container for copying or transferring objects of type <T>.
#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.
A face is a list of labels corresponding to mesh vertices.
const double e
Elementary charge.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Encapsulation of data needed to search for faces.
Cell-face mesh analysis engine.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const Point & rawPoint() const
Return point with no checking.
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Does triangle intersect bounding box.
volumeType getVolumeType(const indexedOctree< treeDataFace > &, const point &) const
Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual const pointField & points() const =0
Return mesh points.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
const vectorField & faceCentres() const
pointField shapePoints() const
Get representative point cloud for all shapes inside.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
treeDataFace(const bool cacheBb, const primitiveMesh &, const labelUList &)
Construct from mesh and subset of faces.
scalar distance() const
Return distance to hit.
direction posBits(const point &) const
Position of point relative to bounding box.
vectorField pointField
pointField is a vectorField.
label start() const
Return start label of this patch in the polyMesh face list.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
errorManip< error > abort(error &err)
const polyMesh & mesh() const
Return the mesh reference.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
label start() const
Return start vertex label.
bool containsAny(const UList< point > &) const
Contains any of the points? (inside or on edge)
defineTypeNameAndDebug(combustionModel, 0)
pointHit nearestPoint(const point &p, const pointField &) const
Return nearest point to face.
bool overlaps(const label index, const treeBoundBox &sampleBb) const
Does (bb of) shape at index overlap bb.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
vector point
Point is a vector.
Non-pointer based hierarchical recursive searching.
const Point & hitPoint() const
Return hit point.
findNearestOp(const indexedOctree< treeDataFace > &tree)
label index() const
Return index.
virtual const faceList & faces() const =0
Return faces.
Standard boundBox + extra functionality for use in octree.
dimensioned< scalar > mag(const dimensioned< Type > &)
bool hit() const
Is there a hit.
A patch is a list of labels that address the faces in the global face list.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
findIntersectOp(const indexedOctree< treeDataFace > &tree)
const labelList & faceLabels() const