boundSphereTemplates.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) 2022-2024 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 "boundSphere.H"
27 #include "labelPair.H"
28 #include "triPointRef.H"
29 #include "tetPointRef.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
34 {
35  return sphere.second() >= 0;
36 }
37 
38 
39 template<class PointField>
42 (
43  const PointField& ps,
44  const FixedList<label, 4>& pis,
45  const label nPs
46 )
47 {
48  switch (nPs)
49  {
50  case 0:
51  return Tuple2<point, scalar>(point::uniform(NaN), -vGreat);
52 
53  case 1:
54  return Tuple2<point, scalar>(ps[pis[0]], 0);
55 
56  case 2:
57  // Return a mid point-centred sphere
58  return Tuple2<point, scalar>
59  (
60  (ps[pis[0]] + ps[pis[1]])/2,
61  mag(ps[pis[0]] - ps[pis[1]])/2
62  );
63 
64  case 3:
65  {
66  // Return the sphere that intersects the triangle
67  return
69  (
70  ps[pis[0]],
71  ps[pis[1]],
72  ps[pis[2]]
73  ).circumCircle();
74  }
75 
76  case 4:
77  {
78  // Return the sphere that intersects the tetrahedron
79  return
81  (
82  ps[pis[0]],
83  ps[pis[1]],
84  ps[pis[2]],
85  ps[pis[3]]
86  ).circumSphere();
87  }
88 
89  default:
91  << "Cannot compute the intersect bounding sphere of more than "
92  << "four points" << exit(FatalError);
93  return Tuple2<point, scalar>(point::uniform(NaN), NaN);
94  }
95 }
96 
97 
98 template<class PointField>
101 (
102  const PointField& ps,
103  const FixedList<label, 4>& pis,
104  const label nPs
105 )
106 {
107  static const scalar tol = 1 + 2*sqrt(small*rootSmall);
108 
109  // Search for spheres that intersect sub-sets of the points that bound all
110  // the other points
111  switch (nPs)
112  {
113  case 0:
114  case 1:
115  case 2:
116  break;
117 
118  case 3:
119  {
120  // Test the spheres that intersect the edges
121  for (label i = 0; i < 3; ++ i)
122  {
123  const point& p0 = ps[pis[i]];
124  const point& p1 = ps[pis[(i + 1) % 3]];
125  const point& pOpp = ps[pis[(i + 2) % 3]];
126 
127  const point c = (p0 + p1)/2;
128  const scalar rSqr = magSqr(p0 - p1)/4;
129 
130  if (magSqr(pOpp - c) <= tol*rSqr)
131  {
132  return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
133  }
134  }
135 
136  break;
137  }
138 
139  case 4:
140  {
141  // Test the spheres that intersect the edges
142  static const FixedList<labelPair, 6> p01s =
143  {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
144  static const FixedList<labelPair, 6> pOpp01s =
145  {{2, 3}, {1, 3}, {1, 2}, {0, 3}, {0, 2}, {0, 1}};
146  for (label i = 0; i < 6; ++ i)
147  {
148  const point& p0 = ps[pis[p01s[i].first()]];
149  const point& p1 = ps[pis[p01s[i].second()]];
150  const point& pOpp0 = ps[pis[pOpp01s[i].first()]];
151  const point& pOpp1 = ps[pis[pOpp01s[i].second()]];
152 
153  const point c = (p0 + p1)/2;
154  const scalar rSqr = magSqr(p0 - p1)/4;
155 
156  if
157  (
158  magSqr(pOpp0 - c) <= tol*rSqr
159  && magSqr(pOpp1 - c) <= tol*rSqr
160  )
161  {
162  return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
163  }
164  }
165 
166  // Test the spheres that intersect the triangles
167  point minC = point::uniform(NaN);
168  scalar minR = vGreat;
169  for (label i = 0; i < 4; ++ i)
170  {
171  const triPointRef tri
172  (
173  ps[pis[i]],
174  ps[pis[(i + 1) % 4]],
175  ps[pis[(i + 2) % 4]]
176  );
177  const point& pOpp = ps[pis[(i + 3) % 4]];
178 
179  const Tuple2<point, scalar> circ = tri.circumCircle();
180  const point& c = circ.first();
181  const scalar r = circ.second();
182 
183  if (magSqr(pOpp - c) <= tol*sqr(r) && r < minR)
184  {
185  minC = c;
186  minR = r;
187  }
188  }
189  if (minR != vGreat)
190  {
191  return Tuple2<point, scalar>(minC, minR);
192  }
193 
194  break;
195  }
196 
197  default:
199  << "Cannot compute the trivial bounding sphere of more than "
200  << "four points" << exit(FatalError);
201  return Tuple2<point, scalar>(point::uniform(NaN), NaN);
202  }
203 
204  // Return the sphere that intersects the points
205  return intersectBoundSphere(ps, pis, nPs);
206 }
207 
208 
209 template<class PointField>
211 (
212  const PointField& ps,
213  List<label>& pis,
214  const label nPs,
215  FixedList<label, 4>& boundaryPis,
216  const label nBoundaryPs
217 )
218 {
219  static const scalar tol = 1 + 2*sqrt(small*rootSmall);
220 
221  Tuple2<point, scalar> sphere;
222 
223  if (nBoundaryPs != 0)
224  {
225  sphere = intersectBoundSphere(ps, boundaryPis, nBoundaryPs);
226  }
227 
228  if (nBoundaryPs == 4)
229  {
230  return sphere;
231  }
232 
233  for (label i = 0; i < nPs; ++ i)
234  {
235  if
236  (
237  (nBoundaryPs == 0 && i == 0)
238  || (magSqr(ps[pis[i]] - sphere.first()) > tol*sqr(sphere.second()))
239  )
240  {
241  boundaryPis[nBoundaryPs] = pis[i];
242 
243  sphere = weizlBoundSphere(ps, pis, i, boundaryPis, nBoundaryPs + 1);
244 
245  // Move the limiting point to the start of the list so that the
246  // sphere grows as quickly as possible in the recursive calls
247  Swap(pis[0], pis[i]);
248  }
249  }
250 
251  return sphere;
252 }
253 
254 
255 template<class PointField>
257 Foam::boundSphere(const PointField& ps, randomGenerator& rndGen)
258 {
259  if (ps.size() <= 4)
260  {
261  static const FixedList<label, 4> pis({0, 1, 2, 3});
262  return trivialBoundSphere(ps, pis, ps.size());
263  }
264  else
265  {
266  labelList pis(identityMap(ps.size()));
267  rndGen.permute(pis);
268  FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
269  return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
270  }
271 }
272 
273 
274 template<class PointField>
276 Foam::boundSphere(const PointField& ps)
277 {
278  if (ps.size() <= 4)
279  {
280  static const FixedList<label, 4> pis({0, 1, 2, 3});
281  return trivialBoundSphere(ps, pis, ps.size());
282  }
283  else
284  {
285  labelList pis(identityMap(ps.size()));
286  randomGenerator(0).permute(pis);
287  FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
288  return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
289  }
290 }
291 
292 
293 // ************************************************************************* //
294 
Functions for constructing bounding spheres of lists of points.
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:66
const Type2 & second() const
Return second.
Definition: Tuple2.H:131
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const dimensionedScalar c
Speed of light in a vacuum.
bool isValidBoundSphere(const Tuple2< point, scalar > &sphere)
Return whether or not the given sphere is valid.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
vector point
Point is a vector.
Definition: point.H:41
Tuple2< point, scalar > intersectBoundSphere(const PointField &ps, const FixedList< label, 4 > &pis, const label nPs)
Compute a sphere of four points or less where every point intersects the.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
Tuple2< point, scalar > trivialBoundSphere(const PointField &ps, const FixedList< label, 4 > &pis, const label nPs)
Compute a bounding sphere of four points or less.
GeometricField< Type, pointPatchField, pointMesh > PointField
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
Tuple2< point, scalar > weizlBoundSphere(const PointField &ps, List< label > &pis, const label nPs, FixedList< label, 4 > &boundaryPis, const label nBoundaryPs)
Compute a bounding sphere for an arbitrary number of points recursively.
error FatalError
void Swap(T &a, T &b)
Definition: Swap.H:43
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Tuple2< point, scalar > boundSphere(const PointField &ps, randomGenerator &rndGen)
Compute a bounding sphere for an arbitrary number of points, and given an.
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
randomGenerator rndGen(653213)