searchableExtrudedCircle.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) 2016-2021 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 
28 #include "Time.H"
29 #include "edgeMesh.H"
30 #include "indexedOctree.H"
31 #include "treeDataEdge.H"
33 #include "quaternion.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(searchableExtrudedCircle, 0);
41  (
42  searchableSurface,
43  searchableExtrudedCircle,
44  dict
45  );
46 }
47 
48 
49 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
50 
52 (
53  const IOobject& io,
54  const dictionary& dict
55 )
56 :
58  eMeshPtr_
59  (
61  (
62  IOobject
63  (
64  dict.lookup("file"),
65  io.time().constant(),
67  io.time(),
70  false
71  ).objectPath(global())
72  )
73  ),
74  radius_(dict.lookup<scalar>("radius"))
75 {
76  const edgeMesh& eMesh = eMeshPtr_();
77 
78  const pointField& points = eMesh.points();
79  const edgeList& edges = eMesh.edges();
80  bounds() = boundBox(points, false);
81 
82  vector halfSpan(0.5*bounds().span());
83  point ctr(bounds().midpoint());
84 
85  bounds().min() = ctr - mag(halfSpan)*vector(1, 1, 1);
86  bounds().max() = ctr + mag(halfSpan)*vector(1, 1, 1);
87 
88  // Calculate bb of all points
89  treeBoundBox bb(bounds());
90 
91  // Slightly extended bb. Slightly off-centred just so on symmetric
92  // geometry there are less face/edge aligned items.
93  bb.min() -= point(rootVSmall, rootVSmall, rootVSmall);
94  bb.max() += point(rootVSmall, rootVSmall, rootVSmall);
95 
96  edgeTree_.reset
97  (
99  (
101  (
102  false, // do not cache bb
103  edges,
104  points,
105  identity(edges.size())
106  ),
107  bb, // overall search domain
108  8, // maxLevel
109  10, // leafsize
110  3.0 // duplicity
111  )
112  );
113 }
114 
115 
116 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
117 
119 {}
120 
121 
122 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
123 
125 {
126  if (regions_.empty())
127  {
128  regions_.setSize(1);
129  regions_[0] = "region0";
130  }
131  return regions_;
132 }
133 
134 
136 {
137  return eMeshPtr_().points().size();
138 }
139 
140 
142 {
143  return eMeshPtr_().points();
144 }
145 
146 
148 (
149  pointField& centres,
150  scalarField& radiusSqr
151 ) const
152 {
153  centres = eMeshPtr_().points();
154  radiusSqr.setSize(centres.size());
155  radiusSqr = Foam::sqr(radius_);
156  // Add a bit to make sure all points are tested inside
157  radiusSqr += Foam::sqr(small);
158 }
159 
160 
162 (
163  const pointField& samples,
164  const scalarField& nearestDistSqr,
165  List<pointIndexHit>& info
166 ) const
167 {
168  const indexedOctree<treeDataEdge>& tree = edgeTree_();
169 
170  info.setSize(samples.size());
171 
172  forAll(samples, i)
173  {
174  info[i] = tree.findNearest(samples[i], nearestDistSqr[i]);
175 
176  if (info[i].hit())
177  {
178  vector d(samples[i]-info[i].hitPoint());
179  info[i].setPoint(info[i].hitPoint() + d/mag(d)*radius_);
180  }
181  }
182 }
183 
184 
186 (
187  const point& start,
188  const point& end,
189  const scalarField& rawLambdas,
190  const scalarField& nearestDistSqr,
191  List<pointIndexHit>& info
192 ) const
193 {
194  const edgeMesh& mesh = eMeshPtr_();
195  const indexedOctree<treeDataEdge>& tree = edgeTree_();
196  const edgeList& edges = mesh.edges();
197  const pointField& points = mesh.points();
198  const labelListList& pointEdges = mesh.pointEdges();
199 
200  const scalar maxDistSqr(Foam::magSqr(bounds().span()));
201 
202  // Normalise lambdas
203  const scalarField lambdas
204  (
205  (rawLambdas-rawLambdas[0])
206  /(rawLambdas.last()-rawLambdas[0])
207  );
208 
209 
210  // Nearest point on curve and local axis direction
211  pointField curvePoints(lambdas.size());
212  vectorField axialVecs(lambdas.size());
213 
214  const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr);
215  curvePoints[0] = startInfo.hitPoint();
216  axialVecs[0] = edges[startInfo.index()].vec(points);
217 
218  const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr);
219  curvePoints.last() = endInfo.hitPoint();
220  axialVecs.last() = edges[endInfo.index()].vec(points);
221 
222 
223 
224  scalarField curveLambdas(points.size(), -1.0);
225 
226  {
227  scalar endDistance = -1.0;
228 
229  // Determine edge lengths walking from start to end.
230 
231  const point& start = curvePoints[0];
232  const point& end = curvePoints.last();
233 
234  label edgei = startInfo.index();
235  const edge& startE = edges[edgei];
236 
237  label pointi = startE[0];
238  if ((startE.vec(points)&(end-start)) > 0)
239  {
240  pointi = startE[1];
241  }
242 
243  curveLambdas[pointi] = mag(points[pointi]-curvePoints[0]);
244  label otherPointi = startE.otherVertex(pointi);
245  curveLambdas[otherPointi] = -mag(points[otherPointi]-curvePoints[0]);
246 
247  // Pout<< "for point:" << points[pointi] << " have distance "
248  // << curveLambdas[pointi] << endl;
249 
250 
251  while (true)
252  {
253  const labelList& pEdges = pointEdges[pointi];
254  if (pEdges.size() == 1)
255  {
256  break;
257  }
258  else if (pEdges.size() != 2)
259  {
260  FatalErrorInFunction << "Curve " << name()
261  << " is not a single path. This is not supported"
262  << exit(FatalError);
263  break;
264  }
265 
266  label oldEdgei = edgei;
267  if (pEdges[0] == oldEdgei)
268  {
269  edgei = pEdges[1];
270  }
271  else
272  {
273  edgei = pEdges[0];
274  }
275 
276  if (edgei == endInfo.index())
277  {
278  endDistance = curveLambdas[pointi] + mag(end-points[pointi]);
279 
280  // Pout<< "Found end edge:" << edges[edgei].centre(points)
281  // << " endPt:" << end
282  // << " point before:" << points[pointi]
283  // << " accumulated length:" << endDistance << endl;
284  }
285 
286 
287  label oldPointi = pointi;
288  pointi = edges[edgei].otherVertex(oldPointi);
289 
290  if (curveLambdas[pointi] >= 0)
291  {
292  break;
293  }
294 
295  curveLambdas[pointi] =
296  curveLambdas[oldPointi] + edges[edgei].mag(points);
297  }
298 
299  // Normalise curveLambdas
300  forAll(curveLambdas, i)
301  {
302  if (curveLambdas[i] >= 0)
303  {
304  curveLambdas[i] /= endDistance;
305  }
306  }
307  }
308 
309 
310 
311  // Interpolation engine
312  linearInterpolationWeights interpolator(curveLambdas);
313 
314  // Find wanted location along curve
315  labelList indices;
316  scalarField weights;
317  for (label i = 1; i < curvePoints.size()-1; i++)
318  {
319  interpolator.valueWeights(lambdas[i], indices, weights);
320 
321  if (indices.size() == 1)
322  {
323  // On outside of curve. Choose one of the connected edges.
324  label pointi = indices[0];
325  const point& p0 = points[pointi];
326  label edge0 = pointEdges[pointi][0];
327  const edge& e0 = edges[edge0];
328  axialVecs[i] = e0.vec(points);
329  curvePoints[i] = weights[0]*p0;
330  }
331  else if (indices.size() == 2)
332  {
333  const point& p0 = points[indices[0]];
334  const point& p1 = points[indices[1]];
335  axialVecs[i] = p1-p0;
336  curvePoints[i] = weights[0]*p0+weights[1]*p1;
337  }
338  }
339  axialVecs /= mag(axialVecs);
340 
341 
342  // Now we have the lambdas, curvePoints and axialVecs.
343 
344 
345 
346  info.setSize(lambdas.size());
347  info = pointIndexHit();
348 
349  // Given the current lambdas interpolate radial direction in between
350  // endpoints (all projected onto the starting coordinate system)
351  quaternion qStart;
352  vector radialStart;
353  {
354  radialStart = start-curvePoints[0];
355  radialStart -= (radialStart&axialVecs[0])*axialVecs[0];
356  radialStart /= mag(radialStart);
357  qStart = quaternion(radialStart, 0.0);
358 
359  info[0] = pointIndexHit(true, start, 0);
360  }
361 
362  quaternion qProjectedEnd;
363  {
364  vector radialEnd(end-curvePoints.last());
365  radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last();
366  radialEnd /= mag(radialEnd);
367 
368  vector projectedEnd = radialEnd;
369  projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0];
370  projectedEnd /= mag(projectedEnd);
371  qProjectedEnd = quaternion(projectedEnd, 0.0);
372 
373  info.last() = pointIndexHit(true, end, 0);
374  }
375 
376  for (label i = 1; i < lambdas.size()-1; i++)
377  {
378  quaternion q(slerp(qStart, qProjectedEnd, lambdas[i]));
379  vector radialDir(q.transform(radialStart));
380 
381  radialDir -= (radialDir&axialVecs[i])*axialVecs.last();
382  radialDir /= mag(radialDir);
383 
384  info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0);
385  }
386 }
387 
388 
390 (
391  const List<pointIndexHit>& info,
392  labelList& region
393 ) const
394 {
395  region.setSize(info.size());
396  region = 0;
397 }
398 
399 
401 (
402  const List<pointIndexHit>& info,
403  vectorField& normal
404 ) const
405 {
406  const edgeMesh& mesh = eMeshPtr_();
407  const indexedOctree<treeDataEdge>& tree = edgeTree_();
408  const edgeList& edges = mesh.edges();
409  const pointField& points = mesh.points();
410 
411  normal.setSize(info.size());
412  normal = Zero;
413 
414  forAll(info, i)
415  {
416  if (info[i].hit())
417  {
418  // Find nearest on curve
419  pointIndexHit curvePt = tree.findNearest
420  (
421  info[i].hitPoint(),
422  Foam::magSqr(bounds().span())
423  );
424 
425  normal[i] = info[i].hitPoint()-curvePt.hitPoint();
426 
427  // Subtract axial direction
428  vector axialVec = edges[curvePt.index()].vec(points);
429  axialVec /= mag(axialVec);
430  normal[i] -= (normal[i]&axialVec)*axialVec;
431  normal[i] /= mag(normal[i]);
432  }
433  }
434 }
435 
436 
437 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual ~searchableExtrudedCircle()
Destructor.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual label size() const
Range of local indices that can be returned.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
static autoPtr< edgeMesh > New(const fileName &, const word &ext)
Select constructed from filename (explicit extension)
Definition: edgeMeshNew.C:31
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Holds data for octree to work on an edges subset.
Definition: treeDataEdge.H:53
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:140
fvMesh & mesh
Macros for easy insertion into run-time selection tables.
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.
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
const Point & hitPoint() const
Return hit point.
const labelListList & pointEdges() const
Return edges.
Definition: edgeMeshI.H:74
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static const word & geometryDir()
Return the geometry directory name.
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
const edgeList & edges() const
Return edges.
Definition: edgeMeshI.H:68
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:60
static const zero Zero
Definition: zero.H:97
dimensioned< scalar > magSqr(const dimensioned< Type > &)
vector vec(const pointField &) const
Return the vector (end - start)
Definition: edgeI.H:175
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
virtual const wordList & regions() const
Names of regions.
const pointField & points() const
Return points.
Definition: edgeMeshI.H:62
defineTypeNameAndDebug(combustionModel, 0)
Points connected by edges.
Definition: edgeMesh.H:69
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
const Time & time() const
Return time.
Definition: IOobject.C:318
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:55
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
searchableExtrudedCircle(const IOobject &io, const dictionary &dict)
Construct from dictionary (used by searchableSurface)
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
A class for managing temporary objects.
Definition: PtrList.H:53
T & last()
Return the last element of the list.
Definition: UListI.H:128
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
label index() const
Return index.
vector transform(const vector &v) const
Rotate the given vector.
Definition: quaternionI.H:296
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864