searchableExtrudedCircle.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) 2016-2017 OpenFOAM Foundation
9  Copyright (C) 2019-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 
31 #include "Time.H"
32 #include "edgeMesh.H"
33 #include "indexedOctree.H"
34 #include "treeDataEdge.H"
36 #include "quaternion.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(searchableExtrudedCircle, 0);
44  (
45  searchableSurface,
46  searchableExtrudedCircle,
47  dict
48  );
50  (
51  searchableSurface,
52  searchableExtrudedCircle,
53  dict,
54  extrudedCircle
55  );
56 }
57 
58 
59 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
60 
61 Foam::searchableExtrudedCircle::searchableExtrudedCircle
62 (
63  const IOobject& io,
64  const dictionary& dict
65 )
66 :
68  eMeshPtr_
69  (
71  (
72  IOobject
73  (
74  dict.get<word>("file"), // name
75  io.time().constant(), // instance
76  "geometry", // local
77  io.time(), // registry
78  IOobject::MUST_READ,
79  IOobject::NO_WRITE,
80  false
81  ).objectPath()
82  )
83  ),
84  radius_(dict.get<scalar>("radius"))
85 {
86  const edgeMesh& eMesh = eMeshPtr_();
87 
88  const pointField& points = eMesh.points();
89  const edgeList& edges = eMesh.edges();
90  bounds() = boundBox(points, false);
91 
92  vector halfSpan(0.5*bounds().span());
93  point ctr(bounds().centre());
94 
95  bounds().min() = ctr - mag(halfSpan) * vector::one;
96  bounds().max() = ctr + mag(halfSpan) * vector::one;
97 
98  // Calculate bb of all points
99  treeBoundBox bb(bounds());
100 
101  // Slightly extended bb. Slightly off-centred just so on symmetric
102  // geometry there are less face/edge aligned items.
103  bb.min() -= point::uniform(ROOTVSMALL);
104  bb.max() += point::uniform(ROOTVSMALL);
105 
106  edgeTree_.reset
107  (
109  (
111  (
112  false, // do not cache bb
113  edges,
114  points,
115  identity(edges.size())
116  ),
117  bb, // overall search domain
118  8, // maxLevel
119  10, // leafsize
120  3.0 // duplicity
121  )
122  );
123 }
124 
125 
126 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
127 
129 {}
130 
131 
132 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
133 
135 {
136  if (regions_.empty())
137  {
138  regions_.setSize(1);
139  regions_[0] = "region0";
140  }
141  return regions_;
142 }
143 
144 
146 {
147  return eMeshPtr_().points().size();
148 }
149 
150 
152 {
153  return eMeshPtr_().points();
154 }
155 
156 
158 (
159  pointField& centres,
160  scalarField& radiusSqr
161 ) const
162 {
163  centres = eMeshPtr_().points();
164  radiusSqr.setSize(centres.size());
165  radiusSqr = Foam::sqr(radius_);
166  // Add a bit to make sure all points are tested inside
167  radiusSqr += Foam::sqr(SMALL);
168 }
169 
170 
172 (
173  const pointField& samples,
174  const scalarField& nearestDistSqr,
175  List<pointIndexHit>& info
176 ) const
177 {
178  const indexedOctree<treeDataEdge>& tree = edgeTree_();
179 
180  info.setSize(samples.size());
181 
182  forAll(samples, i)
183  {
184  const scalar nearestDist = Foam::sqrt(nearestDistSqr[i]);
185  const scalar searchDistSqr = Foam::sqr(nearestDist+radius_);
186 
187  // Find nearest on central edge
188  info[i] = tree.findNearest(samples[i], searchDistSqr);
189 
190  if (info[i].hit())
191  {
192  // Derive distance to nearest surface from distance to nearest edge
193  const vector d(samples[i] - info[i].hitPoint());
194  const scalar s(mag(d));
195 
196  if (s < ROOTVSMALL)
197  {
198  // Point is on edge. TBD.
199  info[i].setMiss();
200  }
201  else
202  {
203  const scalar distToSurface = radius_-s;
204  if (mag(distToSurface) > nearestDist)
205  {
206  info[i].setMiss();
207  }
208  else
209  {
210  info[i].setPoint(info[i].hitPoint() + d/s*radius_);
211  }
212  }
213  }
214  }
215 }
216 
217 
219 (
220  const point& start,
221  const point& end,
222  const scalarField& rawLambdas,
223  const scalarField& nearestDistSqr,
224  List<pointIndexHit>& info
225 ) const
226 {
227  const edgeMesh& mesh = eMeshPtr_();
228  const indexedOctree<treeDataEdge>& tree = edgeTree_();
229  const edgeList& edges = mesh.edges();
230  const pointField& points = mesh.points();
231  const labelListList& pointEdges = mesh.pointEdges();
232 
233  const scalar maxDistSqr(Foam::magSqr(bounds().span()));
234 
235  // Normalise lambdas
236  const scalarField lambdas
237  (
238  (rawLambdas-rawLambdas[0])
239  /(rawLambdas.last()-rawLambdas[0])
240  );
241 
242 
243  // Nearest point on curve and local axis direction
244  pointField curvePoints(lambdas.size());
245  vectorField axialVecs(lambdas.size());
246 
247  const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr);
248  curvePoints[0] = startInfo.hitPoint();
249  axialVecs[0] = edges[startInfo.index()].vec(points);
250 
251  const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr);
252  curvePoints.last() = endInfo.hitPoint();
253  axialVecs.last() = edges[endInfo.index()].vec(points);
254 
255 
256 
257  scalarField curveLambdas(points.size(), -1.0);
258 
259  {
260  scalar endDistance = -1.0;
261 
262  // Determine edge lengths walking from start to end.
263 
264  const point& start = curvePoints[0];
265  const point& end = curvePoints.last();
266 
267  label edgei = startInfo.index();
268  const edge& startE = edges[edgei];
269 
270  label pointi = startE[0];
271  if ((startE.vec(points)&(end-start)) > 0)
272  {
273  pointi = startE[1];
274  }
275 
276  curveLambdas[pointi] = mag(points[pointi]-curvePoints[0]);
277  label otherPointi = startE.otherVertex(pointi);
278  curveLambdas[otherPointi] = -mag(points[otherPointi]-curvePoints[0]);
279 
280  //Pout<< "for point:" << points[pointi] << " have distance "
281  // << curveLambdas[pointi] << endl;
282 
283 
284  while (true)
285  {
286  const labelList& pEdges = pointEdges[pointi];
287  if (pEdges.size() == 1)
288  {
289  break;
290  }
291  else if (pEdges.size() != 2)
292  {
293  FatalErrorInFunction << "Curve " << name()
294  << " is not a single path. This is not supported"
295  << exit(FatalError);
296  break;
297  }
298 
299  label oldEdgei = edgei;
300  if (pEdges[0] == oldEdgei)
301  {
302  edgei = pEdges[1];
303  }
304  else
305  {
306  edgei = pEdges[0];
307  }
308 
309  if (edgei == endInfo.index())
310  {
311  endDistance = curveLambdas[pointi] + mag(end-points[pointi]);
312 
313  //Pout<< "Found end edge:" << edges[edgei].centre(points)
314  // << " endPt:" << end
315  // << " point before:" << points[pointi]
316  // << " accumulated length:" << endDistance << endl;
317  }
318 
319 
320  label oldPointi = pointi;
321  pointi = edges[edgei].otherVertex(oldPointi);
322 
323  if (curveLambdas[pointi] >= 0)
324  {
325  break;
326  }
327 
328  curveLambdas[pointi] =
329  curveLambdas[oldPointi] + edges[edgei].mag(points);
330  }
331 
332  // Normalise curveLambdas
333  forAll(curveLambdas, i)
334  {
335  if (curveLambdas[i] >= 0)
336  {
337  curveLambdas[i] /= endDistance;
338  }
339  }
340  }
341 
342 
343 
344  // Interpolation engine
345  linearInterpolationWeights interpolator(curveLambdas);
346 
347  // Find wanted location along curve
348  labelList indices;
349  scalarField weights;
350  for (label i = 1; i < curvePoints.size()-1; i++)
351  {
352  interpolator.valueWeights(lambdas[i], indices, weights);
353 
354  if (indices.size() == 1)
355  {
356  // On outside of curve. Choose one of the connected edges.
357  label pointi = indices[0];
358  const point& p0 = points[pointi];
359  label edge0 = pointEdges[pointi][0];
360  const edge& e0 = edges[edge0];
361  axialVecs[i] = e0.vec(points);
362  curvePoints[i] = weights[0]*p0;
363  }
364  else if (indices.size() == 2)
365  {
366  const point& p0 = points[indices[0]];
367  const point& p1 = points[indices[1]];
368  axialVecs[i] = p1-p0;
369  curvePoints[i] = weights[0]*p0+weights[1]*p1;
370  }
371  }
372  axialVecs /= mag(axialVecs);
373 
374 
375  // Now we have the lambdas, curvePoints and axialVecs.
376 
377 
378 
379  info.setSize(lambdas.size());
380  info = pointIndexHit();
381 
382  // Given the current lambdas interpolate radial direction inbetween
383  // endpoints (all projected onto the starting coordinate system)
384  quaternion qStart;
385  vector radialStart;
386  {
387  radialStart = start-curvePoints[0];
388  radialStart -= (radialStart&axialVecs[0])*axialVecs[0];
389  radialStart.normalise();
390 
391  qStart = quaternion(radialStart, 0.0);
392 
393  info[0] = pointIndexHit(true, start, 0);
394  }
395 
396  quaternion qProjectedEnd;
397  {
398  vector radialEnd(end-curvePoints.last());
399  radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last();
400  radialEnd.normalise();
401 
402  vector projectedEnd = radialEnd;
403  projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0];
404  projectedEnd.normalise();
405 
406  qProjectedEnd = quaternion(projectedEnd, 0.0);
407 
408  info.last() = pointIndexHit(true, end, 0);
409  }
410 
411  for (label i = 1; i < lambdas.size()-1; i++)
412  {
413  quaternion q(slerp(qStart, qProjectedEnd, lambdas[i]));
414  vector radialDir(q.transform(radialStart));
415 
416  radialDir -= (radialDir & axialVecs[i]) * axialVecs.last();
417  radialDir.normalise();
418 
419  info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0);
420  }
421 }
422 
423 
425 (
426  const List<pointIndexHit>& info,
427  labelList& region
428 ) const
429 {
430  region.setSize(info.size());
431  region = 0;
432 }
433 
434 
436 (
437  const List<pointIndexHit>& info,
438  vectorField& normal
439 ) const
440 {
441  const edgeMesh& mesh = eMeshPtr_();
442  const indexedOctree<treeDataEdge>& tree = edgeTree_();
443  const edgeList& edges = mesh.edges();
444  const pointField& points = mesh.points();
445 
446  normal.setSize(info.size());
447  normal = Zero;
448 
449  forAll(info, i)
450  {
451  if (info[i].hit())
452  {
453  // Find nearest on curve
454  pointIndexHit curvePt = tree.findNearest
455  (
456  info[i].hitPoint(),
457  Foam::magSqr(bounds().span())
458  );
459 
460  normal[i] = info[i].hitPoint()-curvePt.hitPoint();
461 
462  // Subtract axial direction
463  const vector axialVec = edges[curvePt.index()].unitVec(points);
464 
465  normal[i] -= (normal[i] & axialVec) * axialVec;
466  normal[i].normalise();
467  }
468  }
469 }
470 
471 
472 // ************************************************************************* //
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
searchableExtrudedCircle.H
Foam::edgeMesh::points
const pointField & points() const noexcept
Return points.
Definition: edgeMeshI.H:99
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::searchableExtrudedCircle::coordinates
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
Definition: searchableExtrudedCircle.C:151
Foam::searchableExtrudedCircle::findParametricNearest
virtual void findParametricNearest(const point &start, const point &end, const scalarField &lambdas, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Unique to parametric geometry: given points find.
Definition: searchableExtrudedCircle.C:219
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::searchableExtrudedCircle::findNearest
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Definition: searchableExtrudedCircle.C:172
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::slerp
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:78
quaternion.H
Foam::searchableExtrudedCircle::getNormal
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
Definition: searchableExtrudedCircle.C:436
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
indexedOctree.H
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::PointIndexHit::hitPoint
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
Definition: PointIndexHit.H:154
Foam::boundBox::max
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
Foam::primitiveMesh::edges
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Definition: primitiveMeshEdges.C:505
Foam::primitiveMesh::pointEdges
const labelListList & pointEdges() const
Definition: primitiveMeshEdges.C:516
Foam::Vector::normalise
Vector< Cmpt > & normalise()
Normalise the vector by its magnitude.
Definition: VectorI.H:123
Foam::IOobject::time
const Time & time() const
Return Time associated with the objectRegistry.
Definition: IOobject.C:493
Foam::boundBox::min
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::quaternion
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:56
Foam::searchableExtrudedCircle::regions
virtual const wordList & regions() const
Names of regions.
Definition: searchableExtrudedCircle.C:134
Foam::searchableExtrudedCircle::~searchableExtrudedCircle
virtual ~searchableExtrudedCircle()
Destructor.
Definition: searchableExtrudedCircle.C:128
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
linearInterpolationWeights.H
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::Field< vector >
edgeMesh.H
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Definition: searchableSurface.H:69
Foam::addNamedToRunTimeSelectionTable
addNamedToRunTimeSelectionTable(topoSetCellSource, badQualityToCell, word, badQuality)
Foam::indexedOctree
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:50
Foam::searchableExtrudedCircle::getRegion
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
Definition: searchableExtrudedCircle.C:425
Foam::treeDataEdge
Holds data for octree to work on an edges subset.
Definition: treeDataEdge.H:56
samples
scalarField samples(nIntervals, Zero)
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::linearInterpolationWeights::valueWeights
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
Definition: linearInterpolationWeights.C:89
Foam::linearInterpolationWeights
Definition: linearInterpolationWeights.H:50
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
treeDataEdge.H
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::searchableExtrudedCircle::boundingSpheres
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
Definition: searchableExtrudedCircle.C:158
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::edge::otherVertex
label otherVertex(const label pointLabel) const
Given the point label for one vertex, return the other one.
Definition: edgeI.H:188
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Time.H
Foam::edgeMesh::edges
const edgeList & edges() const noexcept
Return edges.
Definition: edgeMeshI.H:105
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::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::edge::vec
vector vec(const UList< point > &pts) const
Return the vector (end - start)
Definition: edgeI.H:417
Foam::Vector< scalar >
Foam::List< edge >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::PointIndexHit::index
label index() const noexcept
Return the hit index.
Definition: PointIndexHit.H:136
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
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::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::TimePaths::constant
const word & constant() const
Return constant name.
Definition: TimePathsI.H:96
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::edgeMesh
Mesh data needed to do the Finite Area discretisation.
Definition: edgeFaMesh.H:53
Foam::searchableExtrudedCircle::size
virtual label size() const
Range of local indices that can be returned.
Definition: searchableExtrudedCircle.C:145
Foam::quaternion::transform
vector transform(const vector &v) const
Rotate the given vector.
Definition: quaternionI.H:327