treeDataFace.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "treeDataFace.H"
27 #include "polyMesh.H"
28 #include "triangleFuncs.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(treeDataFace, 0);
35 
36  scalar treeDataFace::tolSqr = sqr(1e-6);
37 }
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 Foam::treeBoundBox Foam::treeDataFace::calcBb(const label facei) const
43 {
44  const pointField& points = mesh_.points();
45 
46  const face& f = mesh_.faces()[facei];
47 
48  treeBoundBox bb(points[f[0]], points[f[0]]);
49 
50  for (label fp = 1; fp < f.size(); fp++)
51  {
52  const point& p = points[f[fp]];
53 
54  bb.min() = min(bb.min(), p);
55  bb.max() = max(bb.max(), p);
56  }
57  return bb;
58 }
59 
60 
61 void Foam::treeDataFace::update()
62 {
63  forAll(faceLabels_, i)
64  {
65  isTreeFace_.set(faceLabels_[i], 1);
66  }
67 
68  if (cacheBb_)
69  {
70  bbs_.setSize(faceLabels_.size());
71 
72  forAll(faceLabels_, i)
73  {
74  bbs_[i] = calcBb(faceLabels_[i]);
75  }
76  }
77 }
78 
79 
80 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
81 
83 (
84  const bool cacheBb,
85  const primitiveMesh& mesh,
86  const labelUList& faceLabels
87 )
88 :
89  mesh_(mesh),
90  faceLabels_(faceLabels),
91  isTreeFace_(mesh.nFaces(), 0),
92  cacheBb_(cacheBb)
93 {
94  update();
95 }
96 
97 
99 (
100  const bool cacheBb,
101  const primitiveMesh& mesh,
102  labelList&& faceLabels
103 )
104 :
105  mesh_(mesh),
106  faceLabels_(faceLabels),
107  isTreeFace_(mesh.nFaces(), 0),
108  cacheBb_(cacheBb)
109 {
110  update();
111 }
112 
113 
115 (
116  const bool cacheBb,
117  const primitiveMesh& mesh
118 )
119 :
120  mesh_(mesh),
121  faceLabels_(identity(mesh_.nFaces())),
122  isTreeFace_(mesh.nFaces(), 0),
123  cacheBb_(cacheBb)
124 {
125  update();
126 }
127 
128 
130 (
131  const bool cacheBb,
132  const polyPatch& patch
133 )
134 :
135  mesh_(patch.boundaryMesh().mesh()),
136  faceLabels_
137  (
138  identity(patch.size())
139  + patch.start()
140  ),
141  isTreeFace_(mesh_.nFaces(), 0),
142  cacheBb_(cacheBb)
143 {
144  update();
145 }
146 
147 
149 (
150  const indexedOctree<treeDataFace>& tree
151 )
152 :
153  tree_(tree)
154 {}
155 
156 
158 (
159  const indexedOctree<treeDataFace>& tree
160 )
161 :
162  tree_(tree)
163 {}
164 
165 
166 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
167 
169 {
170  pointField cc(faceLabels_.size());
171 
172  forAll(faceLabels_, i)
173  {
174  cc[i] = mesh_.faceCentres()[faceLabels_[i]];
175  }
176 
177  return cc;
178 }
179 
180 
182 (
183  const indexedOctree<treeDataFace>& oc,
184  const point& sample
185 ) const
186 {
187  // Need to determine whether sample is 'inside' or 'outside'
188  // Done by finding nearest face. This gives back a face which is
189  // guaranteed to contain nearest point. This point can be
190  // - in interior of face: compare to face normal
191  // - on edge of face: compare to edge normal
192  // - on point of face: compare to point normal
193  // Unfortunately the octree does not give us back the intersection point
194  // or where on the face it has hit so we have to recreate all that
195  // information.
196 
197 
198  // Find nearest face to sample
199  pointIndexHit info = oc.findNearest(sample, sqr(great));
200 
201  if (info.index() == -1)
202  {
204  << "Could not find " << sample << " in octree."
205  << abort(FatalError);
206  }
207 
208 
209  // Get actual intersection point on face
210  label facei = faceLabels_[info.index()];
211 
212  if (debug & 2)
213  {
214  Pout<< "getSampleType : sample:" << sample
215  << " nearest face:" << facei;
216  }
217 
218  const pointField& points = mesh_.points();
219 
220  // Retest to classify where on face info is. Note: could be improved. We
221  // already have point.
222 
223  const face& f = mesh_.faces()[facei];
224  const vector& area = mesh_.faceAreas()[facei];
225  const point& fc = mesh_.faceCentres()[facei];
226 
227  pointHit curHit = f.nearestPoint(sample, points);
228  const point& curPt = curHit.rawPoint();
229 
230  //
231  // 1] Check whether sample is above face
232  //
233 
234  if (curHit.hit())
235  {
236  // Nearest point inside face. Compare to face normal.
237 
238  if (debug & 2)
239  {
240  Pout<< " -> face hit:" << curPt
241  << " comparing to face normal " << area << endl;
242  }
243  return indexedOctree<treeDataFace>::getSide(area, sample - curPt);
244  }
245 
246  if (debug & 2)
247  {
248  Pout<< " -> face miss:" << curPt;
249  }
250 
251  //
252  // 2] Check whether intersection is on one of the face vertices or
253  // face centre
254  //
255 
256  const scalar typDimSqr = mag(area) + vSmall;
257 
258  forAll(f, fp)
259  {
260  if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
261  {
262  // Face intersection point equals face vertex fp
263 
264  // Calculate point normal (wrong: uses face normals instead of
265  // triangle normals)
266  const labelList& pFaces = mesh_.pointFaces()[f[fp]];
267 
268  vector pointNormal(Zero);
269 
270  forAll(pFaces, i)
271  {
272  if (isTreeFace_.get(pFaces[i]) == 1)
273  {
274  vector n = mesh_.faceAreas()[pFaces[i]];
275  n /= mag(n) + vSmall;
276 
277  pointNormal += n;
278  }
279  }
280 
281  if (debug & 2)
282  {
283  Pout<< " -> face point hit :" << points[f[fp]]
284  << " point normal:" << pointNormal
285  << " distance:"
286  << magSqr(points[f[fp]] - curPt)/typDimSqr << endl;
287  }
289  (
290  pointNormal,
291  sample - curPt
292  );
293  }
294  }
295  if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
296  {
297  // Face intersection point equals face centre. Normal at face centre
298  // is already average of face normals
299 
300  if (debug & 2)
301  {
302  Pout<< " -> centre hit:" << fc
303  << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
304  }
305 
306  return indexedOctree<treeDataFace>::getSide(area, sample - curPt);
307  }
308 
309 
310 
311  //
312  // 3] Get the 'real' edge the face intersection is on
313  //
314 
315  const labelList& myEdges = mesh_.faceEdges()[facei];
316 
317  forAll(myEdges, myEdgeI)
318  {
319  const edge& e = mesh_.edges()[myEdges[myEdgeI]];
320 
321  pointHit edgeHit =
323  (
324  points[e.start()],
325  points[e.end()]
326  ).nearestDist(sample);
327 
328 
329  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
330  {
331  // Face intersection point lies on edge e
332 
333  // Calculate edge normal (wrong: uses face normals instead of
334  // triangle normals)
335  const labelList& eFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
336 
337  vector edgeNormal(Zero);
338 
339  forAll(eFaces, i)
340  {
341  if (isTreeFace_.get(eFaces[i]) == 1)
342  {
343  vector n = mesh_.faceAreas()[eFaces[i]];
344  n /= mag(n) + vSmall;
345 
346  edgeNormal += n;
347  }
348  }
349 
350  if (debug & 2)
351  {
352  Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
353  << " comparing to edge normal:" << edgeNormal
354  << endl;
355  }
356 
357  // Found face intersection point on this edge. Compare to edge
358  // normal
360  (
361  edgeNormal,
362  sample - curPt
363  );
364  }
365  }
366 
367 
368  //
369  // 4] Get the internal edge the face intersection is on
370  //
371 
372  forAll(f, fp)
373  {
375  (
376  points[f[fp]],
377  fc
378  ).nearestDist(sample);
379 
380  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
381  {
382  // Face intersection point lies on edge between two face triangles
383 
384  // Calculate edge normal as average of the two triangle normals
385  vector e = points[f[fp]] - fc;
386  vector ePrev = points[f[f.rcIndex(fp)]] - fc;
387  vector eNext = points[f[f.fcIndex(fp)]] - fc;
388 
389  vector nLeft = ePrev ^ e;
390  nLeft /= mag(nLeft) + vSmall;
391 
392  vector nRight = e ^ eNext;
393  nRight /= mag(nRight) + vSmall;
394 
395  if (debug & 2)
396  {
397  Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
398  << " comparing to edge normal "
399  << 0.5*(nLeft + nRight)
400  << endl;
401  }
402 
403  // Found face intersection point on this edge. Compare to edge
404  // normal
406  (
407  0.5*(nLeft + nRight),
408  sample - curPt
409  );
410  }
411  }
412 
413  if (debug & 2)
414  {
415  Pout<< "Did not find sample " << sample
416  << " anywhere related to nearest face " << facei << endl
417  << "Face:";
418 
419  forAll(f, fp)
420  {
421  Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
422  << endl;
423  }
424  }
425 
426  // Can't determine status of sample with respect to nearest face.
427  // Either
428  // - tolerances are wrong. (if e.g. face has zero area)
429  // - or (more likely) surface is not closed.
430 
431  return volumeType::unknown;
432 }
433 
434 
435 // Check if any point on shape is inside cubeBb.
437 (
438  const label index,
439  const treeBoundBox& cubeBb
440 ) const
441 {
442  // 1. Quick rejection: bb does not intersect face bb at all
443  if (cacheBb_)
444  {
445  if (!cubeBb.overlaps(bbs_[index]))
446  {
447  return false;
448  }
449  }
450  else
451  {
452  if (!cubeBb.overlaps(calcBb(faceLabels_[index])))
453  {
454  return false;
455  }
456  }
457 
458  const pointField& points = mesh_.points();
459 
460 
461  // 2. Check if one or more face points inside
462  label facei = faceLabels_[index];
463 
464  const face& f = mesh_.faces()[facei];
465  if (cubeBb.containsAny(points, f))
466  {
467  return true;
468  }
469 
470  // 3. Difficult case: all points are outside but connecting edges might
471  // go through cube. Use triangle-bounding box intersection.
472  const point& fc = mesh_.faceCentres()[facei];
473 
474  forAll(f, fp)
475  {
476  bool triIntersects = triangleFuncs::intersectBb
477  (
478  points[f[fp]],
479  points[f[f.fcIndex(fp)]],
480  fc,
481  cubeBb
482  );
483 
484  if (triIntersects)
485  {
486  return true;
487  }
488  }
489  return false;
490 }
491 
492 
493 void Foam::treeDataFace::findNearestOp::operator()
494 (
495  const labelUList& indices,
496  const point& sample,
497 
498  scalar& nearestDistSqr,
499  label& minIndex,
500  point& nearestPoint
501 ) const
502 {
503  const treeDataFace& shape = tree_.shapes();
504 
505  forAll(indices, i)
506  {
507  const label index = indices[i];
508 
509  const face& f = shape.mesh().faces()[shape.faceLabels()[index]];
510 
511  pointHit nearHit = f.nearestPoint(sample, shape.mesh().points());
512  scalar distSqr = sqr(nearHit.distance());
513 
514  if (distSqr < nearestDistSqr)
515  {
516  nearestDistSqr = distSqr;
517  minIndex = index;
518  nearestPoint = nearHit.rawPoint();
519  }
520  }
521 }
522 
523 
524 void Foam::treeDataFace::findNearestOp::operator()
525 (
526  const labelUList& indices,
527  const linePointRef& ln,
528 
529  treeBoundBox& tightest,
530  label& minIndex,
531  point& linePoint,
532  point& nearestPoint
533 ) const
534 {
536 }
537 
538 
539 bool Foam::treeDataFace::findIntersectOp::operator()
540 (
541  const label index,
542  const point& start,
543  const point& end,
544  point& intersectionPoint
545 ) const
546 {
547  const treeDataFace& shape = tree_.shapes();
548 
549  // Do quick rejection test
550  if (shape.cacheBb_)
551  {
552  const treeBoundBox& faceBb = shape.bbs_[index];
553 
554  if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
555  {
556  // start and end in same block outside of faceBb.
557  return false;
558  }
559  }
560 
561  const label facei = shape.faceLabels_[index];
562 
563  const vector dir(end - start);
564 
565  pointHit inter = shape.mesh_.faces()[facei].intersection
566  (
567  start,
568  dir,
569  shape.mesh_.faceCentres()[facei],
570  shape.mesh_.points(),
572  );
573 
574  if (inter.hit() && inter.distance() <= 1)
575  {
576  // Note: no extra test on whether intersection is in front of us
577  // since using half_ray
578  intersectionPoint = inter.hitPoint();
579  return true;
580  }
581  else
582  {
583  return false;
584  }
585 }
586 
587 
588 // ************************************************************************* //
const primitiveMesh & mesh() const
Definition: treeDataFace.H:184
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A line primitive.
Definition: line.H:56
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Encapsulation of data needed to search for faces.
Definition: treeDataFace.H:58
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:297
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Does triangle intersect bounding box.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
pointField shapePoints() const
Get representative point cloud for all shapes inside.
Definition: treeDataFace.C:168
virtual const pointField & points() const =0
Return mesh points.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
bool overlaps(const label index, const treeBoundBox &sampleBb) const
Does (bb of) shape at index overlap bb.
Definition: treeDataFace.C:437
direction posBits(const point &) const
Position of point relative to bounding box.
Definition: treeBoundBox.C:465
static const layerAndWeight max
static const Form min
Definition: VectorSpace.H:116
treeDataFace(const bool cacheBb, const primitiveMesh &, const labelUList &)
Construct from mesh and subset of faces.
Definition: treeDataFace.C:83
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
const Point & hitPoint() const
Return hit point.
Definition: PointHit.H:126
const polyMesh & mesh() const
Return the mesh reference.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
bool hit() const
Is there a hit.
Definition: PointHit.H:120
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
Definition: POSIX.C:908
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, &mergedCyclicPolyPatch::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]
Definition: readKivaGrid.H:235
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
defineTypeNameAndDebug(combustionModel, 0)
const vectorField & faceCentres() const
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
pointHit nearestPoint(const point &p, const pointField &) const
Return nearest point to face.
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
findNearestOp(const indexedOctree< treeDataFace > &tree)
Definition: treeDataFace.C:149
virtual const faceList & faces() const =0
Return faces.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
volScalarField & p
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
volumeType getVolumeType(const indexedOctree< treeDataFace > &, const point &) const
Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
Definition: treeDataFace.C:182
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:353
findIntersectOp(const indexedOctree< treeDataFace > &tree)
Definition: treeDataFace.C:158
label index() const
Return index.
label start() const
Return start vertex label.
Definition: edgeI.H:81
const labelList & faceLabels() const
Definition: treeDataFace.H:179
Namespace for OpenFOAM.
bool containsAny(const UList< point > &) const
Contains any of the points? (inside or on edge)
Definition: boundBox.C:261