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