ESyS-Particle  4.0.1
IntersectionVolCalculator.hpp
1 
2 // //
3 // Copyright (c) 2003-2011 by The University of Queensland //
4 // Earth Systems Science Computational Centre (ESSCC) //
5 // http://www.uq.edu.au/esscc //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 
15 #include <stdexcept>
16 #include <sstream>
17 
18 namespace esys
19 {
20  namespace lsm
21  {
22  namespace impl
23  {
24  inline double square(double val)
25  {
26  return val*val;
27  }
28 
29  template <int tmplDim, typename TmplVec>
30  DimBasicBox<tmplDim,TmplVec>::DimBasicBox(const Vec &minPt, const Vec &maxPt)
31  : m_minPt(minPt),
32  m_maxPt(maxPt)
33  {
34  }
35 
36  template <int tmplDim, typename TmplVec>
37  const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMinPt() const
38  {
39  return m_minPt;
40  }
41 
42  template <int tmplDim, typename TmplVec>
43  const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMaxPt() const
44  {
45  return m_maxPt;
46  }
47 
48  template <int tmplDim, typename TmplVec>
49  double DimBasicBox<tmplDim,TmplVec>::getVolume() const
50  {
51  double product = 1.0;
52  for (int i = 0; i < tmplDim; i++)
53  {
54  product *= (this->getMaxPt()[i] - this->getMinPt()[i]);
55  }
56  return product;
57  }
58 
59  template <int tmplDim, typename TmplVec>
60  template <typename TmplSphere>
61  bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const TmplSphere &sphere) const
62  {
63  double distSqrd = 0.0;
64  for (int i = 0; i < tmplDim; i++)
65  {
66  if (sphere.getCentre()[i] < this->getMinPt()[i])
67  {
68  distSqrd += square(sphere.getCentre()[i] - this->getMinPt()[i]);
69  }
70  else if (sphere.getCentre()[i] > this->getMaxPt()[i])
71  {
72  distSqrd += square(sphere.getCentre()[i] - this->getMaxPt()[i]);
73  }
74  }
75  return (distSqrd <= square(sphere.getRadius()));
76  }
77 
78  template <int tmplDim, typename TmplVec>
79  bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const
80  {
81  for (int i = 0; (i < tmplDim); i++)
82  {
83  if ((this->getMinPt()[i] > pt[i]) || (pt[i] > this->getMaxPt()[i]))
84  {
85  return false;
86  }
87  }
88  return true;
89  }
90 
91  template <int tmplDim, typename TmplVec>
92  template <typename TmplSphere>
93  bool DimBasicBox<tmplDim,TmplVec>::contains(const TmplSphere &sphere) const
94  {
95  for (int i = 0; (i < tmplDim); i++)
96  {
97  Vec pt = Vec(0.0);
98  pt[i] = sphere.getRadius();
99  if
100  (
101  !(intersectsWith(sphere.getCentre() + pt))
102  ||
103  !(intersectsWith(sphere.getCentre() - pt))
104  )
105  {
106  return false;
107  }
108  }
109  return true;
110  }
111 
112  template <int tmplDim, typename TmplVec>
113  double DimPlane<tmplDim,TmplVec>::norm(const Vec &pt)
114  {
115  double sum = 0.0;
116  for (int i = 0; i < tmplDim; i++)
117  {
118  sum += pt[i]*pt[i];
119  }
120  return sqrt(sum);
121  }
122 
123  template <int tmplDim, typename TmplVec>
124  double DimPlane<tmplDim,TmplVec>::dot(const Vec &p1, const Vec &p2)
125  {
126  double sum = 0.0;
127  for (int i = 0; i < tmplDim; i++)
128  {
129  sum += p1[i]*p2[i];
130  }
131  return sum;
132  }
133 
134  template <int tmplDim, typename TmplVec>
135  DimPlane<tmplDim,TmplVec>::DimPlane() : m_normal(), m_pt(), m_invNormalNorm(0.0)
136  {
137  }
138 
139  template <int tmplDim, typename TmplVec>
140  DimPlane<tmplDim,TmplVec>::DimPlane(const Vec &normal, const Vec &pt)
141  : m_normal(normal),
142  m_pt(pt),
143  m_invNormalNorm((1.0/norm(normal)))
144  {
145  }
146 
147  template <int tmplDim, typename TmplVec>
148  DimPlane<tmplDim,TmplVec>::DimPlane(const DimPlane &plane)
149  : m_normal(plane.m_normal),
150  m_pt(plane.m_pt),
151  m_invNormalNorm(plane.m_invNormalNorm)
152  {
153  }
154 
155  template <int tmplDim, typename TmplVec>
156  DimPlane<tmplDim,TmplVec> &DimPlane<tmplDim,TmplVec>::operator=(const DimPlane &plane)
157  {
158  m_normal = plane.m_normal;
159  m_pt = plane.m_pt;
160  m_invNormalNorm = plane.m_invNormalNorm;
161  return *this;
162  }
163 
164  template <int tmplDim, typename TmplVec>
165  double DimPlane<tmplDim,TmplVec>::getSignedDistanceTo(const Vec &pt) const
166  {
167  // http://mathworld.wolfram.com/Point-PlaneDistance.html
168  return
169  (
170  (dot(m_normal, pt) - dot(m_normal, m_pt))*m_invNormalNorm
171  );
172  }
173 
174  template <int tmplDim, typename TmplVec>
175  double DimPlane<tmplDim,TmplVec>::getDistanceTo(const Vec &pt) const
176  {
177  return fabs(getSignedDistanceTo(pt));
178  }
179 
180  template <int tmplDim, typename TmplVec>
181  const typename DimPlane<tmplDim,TmplVec>::Vec &DimPlane<tmplDim,TmplVec>::getNormal() const
182  {
183  return m_normal;
184  }
185 
186  template <int tmplDim, typename TmplVec>
187  const double DimBasicSphere<tmplDim,TmplVec>::FOUR_THIRDS_PI = (4.0/3.0)*M_PI;
188 
189  template <int tmplDim, typename TmplVec>
190  const double DimBasicSphere<tmplDim,TmplVec>::ONE_THIRD_PI = M_PI/3.0;
191 
192  template <int tmplDim, typename TmplVec>
193  DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere()
194  : m_centre(),
195  m_radius(0.0)
196  {
197  }
198 
199  template <int tmplDim, typename TmplVec>
200  DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const Vec &centrePt, double radius)
201  : m_centre(centrePt),
202  m_radius(radius)
203  {
204  }
205 
206  template <int tmplDim, typename TmplVec>
207  DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const DimBasicSphere &sphere)
208  : m_centre(sphere.m_centre),
209  m_radius(sphere.m_radius)
210  {
211  }
212 
213  template <int tmplDim, typename TmplVec>
214  DimBasicSphere<tmplDim,TmplVec> &DimBasicSphere<tmplDim,TmplVec>::operator=(const DimBasicSphere &sphere)
215  {
216  m_centre = sphere.m_centre;
217  m_radius = sphere.m_radius;
218  return *this;
219  }
220 
221  template <int tmplDim, typename TmplVec>
222  double DimBasicSphere<tmplDim,TmplVec>::getRadius() const
223  {
224  return m_radius;
225  }
226 
227  template <int tmplDim, typename TmplVec>
228  const typename DimBasicSphere<tmplDim,TmplVec>::Vec &
229  DimBasicSphere<tmplDim,TmplVec>::getCentre() const
230  {
231  return m_centre;
232  }
233 
234  template <int tmplDim, typename TmplVec>
235  double DimBasicSphere<tmplDim,TmplVec>::getVolume() const
236  {
237  return (tmplDim == 2) ? M_PI*getRadius()*getRadius() : FOUR_THIRDS_PI*getRadius()*getRadius()*getRadius();
238  }
239 
240  inline void checkDomain(double r, double x1, double y1, double x2, double y2)
241  {
242  const double rSqrd = r*r;
243  const double x1Sqrd = x1*x1;
244  const double x2Sqrd = x2*x2;
245  const double y1Sqrd = y1*y1;
246  const double y2Sqrd = y2*y2;
247  if
248  (
249  ((rSqrd - x1Sqrd - y1Sqrd) < 0)
250  ||
251  ((rSqrd - x1Sqrd - y2Sqrd) < 0)
252  ||
253  ((rSqrd - x2Sqrd - y1Sqrd) < 0)
254  ||
255  ((rSqrd - x2Sqrd - y2Sqrd) < 0)
256  )
257  {
258  std::stringstream msg;
259  msg
260  << "Invalid rectangular domain for sphere integration, points in domain "
261  << "(" << x1 << "," << y1 << "), (" << x2 << "," << y2 << " lie outside "
262  << "sphere of radius " << r << " centred at the origin.";
263  throw std::runtime_error(msg.str());
264  }
265  }
266 
267  template <int tmplDim, typename TmplVec>
268  double DimBasicSphere<tmplDim,TmplVec>::getVolume(const Vec &minPt, const Vec &maxPt, const int dimX, const int dimY) const
269  {
270  double vol = 0.0;
271  if ((tmplDim == 2) || (tmplDim == 3))
272  {
273  if (minPt[dimX] != maxPt[dimX])
274  {
275  const double x1 = minPt[dimX] - getCentre()[dimX];
276  const double x2 = maxPt[dimX] - getCentre()[dimX];
277  const double r = getRadius();
278 
279  if (tmplDim == 2)
280  {
281  const double rSqrd = r*r;
282  const double x1Sqrd = x1*x1;
283  const double x2Sqrd = x2*x2;
284 
285  vol =
286  (
287  rSqrd*asin(x2/r)
288  +
289  x2*sqrt(rSqrd-x2Sqrd)
290  -
291  rSqrd*asin(x1/r)
292  -
293  x1*sqrt(rSqrd-x1Sqrd)
294  )*0.5;
295  }
296  else if (tmplDim == 3)
297  {
298  if (minPt[dimY] != maxPt[dimY])
299  {
300  const double y1 = minPt[dimY] - getCentre()[dimY];
301  const double y2 = maxPt[dimY] - getCentre()[dimY];
302 
303  checkDomain(r, x1, y1, x2, y2);
304 
305  //
306  // Matlab6/maple generated code:
307  //
308  // syms x y x1 x2 y1 y2 r real
309  // sphereIntegral = int(int('sqrt(r^2-x^2-y^2)', x, x1, x2), y, y1, y2)
310  // maple('readlib(codegen)')
311  // maple('readlib(optimize)')
312  // optimizedIntegral = maple('optimize', sphereIntegral, 'tryhard')
313  // maple('cost', optimizedIntegral)
314  //
315  //43*additions+82*multiplications+6*divisions+22*functions+42*assignments
316 
317  const double t30 = y2*y2; //y2^2;
318  const double t31 = x2*x2; //x2^2;
319  const double t36 = r*r; //r^2;
320  const double t59 = t31-t36;
321  const double t40 = sqrt(-t30-t59); //pow(-t30-t59,1/2);
322  const double t10 = 1.0/t40;
323  const double t32 = x1*x1; //x1^2;
324  const double t54 = t32-t36;
325  const double t42 = sqrt(-t30-t54); //pow(-t30-t54,1/2);
326  const double t14 = 1.0/t42;
327  const double t64 = -atan(t10*x2)+atan(t14*x1);
328  const double t27 = y1*y1; //y1^2;
329  const double t39 = sqrt(-t27-t59); //pow(-t27-t59,1/2);
330  const double t9 = 1.0/t39;
331  const double t41 = sqrt(-t27-t54); //pow(-t27-t54,1/2);
332  const double t12 = 1.0/t41;
333  const double t63 = atan(t12*x1)-atan(t9*x2);
334  const double t62 = -atan(y2*t14)+atan(y1*t12);
335  const double t61 = -atan(y1*t9)+atan(t10*y2);
336  const double t37 = sqrt(t31); //pow(t31,1/2);
337  const double t21 = 1.0/t37;
338  const double t60 = t21*t9;
339  const double t38 = sqrt(t32); //pow(t32,1/2);
340  const double t24 = 1.0/t38;
341  const double t58 = t24*t14;
342  const double t57 = t24*t12;
343  const double t56 = t37*t38;
344  const double t55 = t21*t10;
345  const double t53 = 2.0*x2;
346  const double t52 = 2.0*x1;
347  const double t51 = t42*t56;
348  const double t28 = t27*y1;
349  const double t50 = t28-t36*y1;
350  const double t34 = t30*y2;
351  const double t49 = t34-t36*y2;
352  const double t48 = t41*t51;
353  const double t35 = t36*r;
354  const double t33 = t31*x2;
355  const double t29 = t32*x1;
356  const double t26 = r*y2;
357  const double t25 = y1*r;
358  vol =
359  (-1.0/6.0)
360  *
361  (
362  (-2.0*t33*y1-t50*t53)*t40*t48
363  +
364  (
365  (2.0*t33*y2+t49*t53)*t48
366  +
367  (
368  (2.0*t29*y1+t50*t52)*t51
369  +
370  (
371  (-2.0*t29*y2-t49*t52)*t56
372  +
373  (
374  (
375  atan((t25+t54)*t57)
376  +
377  atan((-t26+t54)*t58)
378  -
379  atan((t26+t54)*t58)
380  -
381  atan((-t25+t54)*t57)
382  )*t35*t37*x1
383  +
384  (
385  (
386  -atan((-t26+t59)*t55)
387  -
388  atan((t25+t59)*t60)
389  +
390  atan((-t25+t59)*t60)
391  +
392  atan((t26+t59)*t55)
393  )*t35*x2
394  +
395  (-t64*t34+t61*t33+t62*t29+t63*t28+3.0*(t64*y2-t63*y1-t61*x2-t62*x1)*t36)*t37
396  )*t38
397  )*t42
398  )*t41
399  )*t40
400  )*t39
401  )*t14*t9*t55*t57;
402  }
403  }
404  }
405  }
406 
407  return vol;
408  }
409 
410  template <int tmplDim, typename TmplVec>
411  bool DimBasicSphere<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const
412  {
413  double distSqrd = 0.0;
414  for (int i = 0; i < tmplDim; i++)
415  {
416  distSqrd += square(getCentre()[i] - pt[i]);
417  }
418  return (distSqrd <= square(getRadius()));
419  }
420 
421  template <int tmplDim, typename TmplVec>
422  double DimBasicSphere<tmplDim,TmplVec>::getSegmentVolume(const Plane &plane) const
423  {
424  double vol = 0.0;
425  if ((tmplDim == 2) || (tmplDim == 3))
426  {
427  const double signedD = plane.getSignedDistanceTo(getCentre());
428  const double d = fabs(signedD);
429  if (d < getRadius())
430  {
431  if (tmplDim == 2)
432  {
433  // http://mathworld.wolfram.com/CircularSegment.html
434  const double rSqrd = getRadius()*getRadius();
435  vol = rSqrd*acos(d/getRadius()) - d*sqrt(rSqrd - d*d);
436  }
437  else if (tmplDim == 3)
438  {
439  // http://mathworld.wolfram.com/SphericalCap.html
440  const double h = getRadius() - d;
441  vol = ONE_THIRD_PI*h*h*(3.0*getRadius()-h);
442  }
443  vol = ((signedD < 0) ? vol : getVolume() - vol);
444  }
445  }
446  return vol;
447  }
448 
449  template <int tmplDim, typename TmplVec>
450  typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
451  IntersectionVolCalculator<tmplDim,TmplVec>::getNormal(int dim)
452  {
453  Vec n = Vec(0.0);
454  n[dim] = 1.0;
455  return n;
456  }
457 
458  template <int tmplDim, typename TmplVec>
459  typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
460  IntersectionVolCalculator<tmplDim,TmplVec>::getNegNormal(int dim)
461  {
462  Vec n = Vec(0.0);
463  n[dim] = -1.0;
464  return n;
465  }
466 
467  template <int tmplDim, typename TmplVec>
468  IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere()
469  : m_sphere(),
470  m_volume(0.0)
471  {
472  }
473 
474  template <int tmplDim, typename TmplVec>
475  IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
476  const BasicSphere &sphere
477  )
478  : m_sphere(sphere),
479  m_volume(sphere.getVolume())
480  {
481  }
482 
483  template <int tmplDim, typename TmplVec>
484  IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
485  const VolumeSphere &sphere
486  )
487  : m_sphere(BasicSphere(sphere.getCentre(), sphere.getRadius())),
488  m_volume(sphere.m_volume)
489  {
490  }
491 
492  template <int tmplDim, typename TmplVec>
493  typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
494  IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::operator=(
495  const VolumeSphere &sphere
496  )
497  {
498  m_sphere = sphere.m_sphere;
499  m_volume = sphere.m_volume;
500  return *this;
501  }
502 
503  template <int tmplDim, typename TmplVec>
504  double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getRadius() const
505  {
506  return m_sphere.getRadius();
507  }
508 
509  template <int tmplDim, typename TmplVec>
510  const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
511  IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getCentre() const
512  {
513  return m_sphere.getCentre();
514  }
515 
516  template <int tmplDim, typename TmplVec>
517  double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume() const
518  {
519  return m_volume;
520  }
521 
522  template <int tmplDim, typename TmplVec>
523  double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume(
524  const Vec &minPt,
525  const Vec &maxPt,
526  const int dimX,
527  const int dimY
528  ) const
529  {
530  return m_sphere.getVolume(minPt, maxPt, dimX, dimY);
531  }
532 
533  template <int tmplDim, typename TmplVec>
534  double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::calcVolume() const
535  {
536  return m_sphere.getVolume();
537  }
538 
539  template <int tmplDim, typename TmplVec>
540  bool IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::intersectsWith(const Vec &pt) const
541  {
542  return m_sphere.intersectsWith(pt);
543  }
544 
545  template <int tmplDim, typename TmplVec>
546  double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getSegmentVolume(const Plane &plane) const
547  {
548  return m_sphere.getSegmentVolume(plane);
549  }
550 
551  template <int tmplDim, typename TmplVec>
552  IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex() : m_pt()
553  {
554  }
555 
556  template <int tmplDim, typename TmplVec>
557  IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vec &pt) : m_pt(pt)
558  {
559  }
560 
561  template <int tmplDim, typename TmplVec>
562  IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vertex &vtx) : m_pt(vtx.m_pt)
563  {
564  }
565 
566  template <int tmplDim, typename TmplVec>
567  typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
568  IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::operator=(const Vertex &vtx)
569  {
570  m_pt = vtx.m_pt;
571  return *this;
572  }
573 
574  template <int tmplDim, typename TmplVec>
575  const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
576  IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::getPoint() const
577  {
578  return m_pt;
579  }
580 
581  template <int tmplDim, typename TmplVec>
582  void IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::setPoint(const Vec &pt)
583  {
584  m_pt = pt;
585  }
586 
587  template <int tmplDim, typename TmplVec>
588  IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
589  const BasicBox &box
590  )
591  : BasicBox(box),
592  m_vertexArray()
593  {
594  createVertices();
595  }
596 
597  template <int tmplDim, typename TmplVec>
598  IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
599  const VertexBox &box
600  )
601  : BasicBox(box)
602  {
603  for (int i = 0; i < getNumVertices(); i++)
604  {
605  m_vertexArray[i] = box.m_vertexArray[i];
606  }
607  }
608 
609  template <int tmplDim, typename TmplVec>
610  typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
611  IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::operator=(
612  const VertexBox &box
613  )
614  {
615  BasicBox::operator=(box);
616  for (int i = 0; i < getNumVertices(); i++)
617  {
618  m_vertexArray[i] = box.m_vertexArray[i];
619  }
620  return *this;
621  }
622 
623  template <int tmplDim, typename TmplVec>
624  void IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::createVertices()
625  {
626  int j = 0;
627  m_vertexArray[j].setPoint(this->getMinPt());
628  int i = 0;
629  for (j++; i < tmplDim; i++, j++)
630  {
631  Vec pt = this->getMinPt();
632  pt[i] = this->getMaxPt()[i];
633  m_vertexArray[j].setPoint(pt);
634  }
635 
636  m_vertexArray[j].setPoint(this->getMaxPt());
637  for (i = 0, j++; i < tmplDim && j < s_numVertices; i++, j++)
638  {
639  Vec pt = this->getMaxPt();
640  pt[i] = this->getMinPt()[i];
641  m_vertexArray[j] = pt;
642  }
643  }
644 
645  template <int tmplDim, typename TmplVec>
646  const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
647  IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getVertex(int i) const
648  {
649  return m_vertexArray[i];
650  }
651 
652  template <int tmplDim, typename TmplVec>
653  int IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getNumVertices()
654  {
655  return s_numVertices;
656  }
657 
658  template <int tmplDim, typename TmplVec>
659  IntersectionVolCalculator<tmplDim,TmplVec>::IntersectionVolCalculator(
660  const BasicBox &box
661  )
662  : m_sphere(BasicSphere(Vec(), 1.0)),
663  m_box(box)
664  {
665  }
666 
667  template <int tmplDim, typename TmplVec>
668  const typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
669  IntersectionVolCalculator<tmplDim,TmplVec>::getSphere() const
670  {
671  return m_sphere;
672  }
673 
674  template <int tmplDim, typename TmplVec>
675  void IntersectionVolCalculator<tmplDim,TmplVec>::setSphere(
676  const BasicSphere &sphere
677  )
678  {
679  m_sphere = sphere;
680  }
681 
682  template <int tmplDim, typename TmplVec>
683  const typename IntersectionVolCalculator<tmplDim,TmplVec>::BasicBox &
684  IntersectionVolCalculator<tmplDim,TmplVec>::getBox() const
685  {
686  return m_box;
687  }
688 
689  template <int tmplDim, typename TmplVec>
690  const typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
691  IntersectionVolCalculator<tmplDim,TmplVec>::getVertexBox() const
692  {
693  return m_box;
694  }
695 
696  template <int tmplDim, typename TmplVec>
697  typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
698  IntersectionVolCalculator<tmplDim,TmplVec>::componentMin(
699  const Vec &p1,
700  const Vec &p2
701  )
702  {
703  Vec m;
704  for (int i = 0; i < tmplDim; i++)
705  {
706  m[i] = min(p1[i], p2[i]);
707  }
708  return m;
709  }
710 
711  template <int tmplDim, typename TmplVec>
712  typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
713  IntersectionVolCalculator<tmplDim,TmplVec>::componentMax(
714  const Vec &p1,
715  const Vec &p2
716  )
717  {
718  Vec m;
719  for (int i = 0; i < tmplDim; i++)
720  {
721  m[i] = max(p1[i], p2[i]);
722  }
723  return m;
724  }
725 
726  template <int tmplDim, typename TmplVec>
727  double IntersectionVolCalculator<tmplDim,TmplVec>::getInsidePointVolume(
728  const Vec &pt
729  ) const
730  {
731  double vol = 0.0;
732  const Vec &centrePt = getSphere().getCentre();
733  const Vec diag = (getSphere().getCentre() - pt)*2.0;
734  const Vec oppCorner = pt + diag;
735  BasicBox box =
736  BasicBox(
737  componentMin(pt, oppCorner),
738  componentMax(pt, oppCorner)
739  );
740  const double boxVol = box.getVolume();
741  const double sphVol = getSphere().getVolume();
742 
743  double s[tmplDim];
744  double v[tmplDim+1];
745  for (int i = 0; i < tmplDim; i++)
746  {
747  s[i] = getSphere().getSegmentVolume(Plane(getNormal((i+1)%tmplDim), box.getMaxPt()));
748  }
749  if (tmplDim == 2)
750  {
751  v[0] = 0.50*(sphVol - 2.0*s[0] - boxVol);
752  v[1] = 0.50*(sphVol - 2.0*s[1] - boxVol);
753  v[2] = 0.25*(sphVol - 2.0*v[0] - 2.0*v[1] - boxVol);
754 
755  if (pt[0] <= centrePt[0])
756  {
757  if (pt[1] <= centrePt[1])
758  {
759  vol = boxVol + v[0] + v[1] + v[2];
760  }
761  else
762  {
763  vol = v[1] + v[2];
764  }
765  }
766  else
767  {
768  if (pt[1] <= centrePt[1])
769  {
770  vol = v[0] + v[2];
771  }
772  else
773  {
774  vol = v[2];
775  }
776  }
777  }
778  else if (tmplDim == 3)
779  {
780  v[0] =
781  0.500*(
782  2.0*getSphere().getVolume(
783  box.getMinPt(),
784  box.getMaxPt(),
785  1,
786  2
787  )
788  -
789  boxVol
790  );
791  v[1] =
792  0.500*(
793  2.0*getSphere().getVolume(
794  box.getMinPt(),
795  box.getMaxPt(),
796  0,
797  2
798  )
799  -
800  boxVol
801  );
802  v[2] =
803  0.500*(
804  2.0*getSphere().getVolume(
805  box.getMinPt(),
806  box.getMaxPt(),
807  0,
808  1
809  )
810  -
811  boxVol
812  );
813 
814  double e[3];
815  e[0] = 0.250*(sphVol - 2.0*s[1] - 2.0*v[0] - 2.0*v[1] - boxVol);
816  e[1] = 0.250*(sphVol - 2.0*s[2] - 2.0*v[1] - 2.0*v[2] - boxVol);
817  e[2] = 0.250*(sphVol - 2.0*s[0] - 2.0*v[0] - 2.0*v[2] - boxVol);
818 
819  v[3] =
820  0.125*(
821  sphVol
822  -
823  2.0*v[0] - 2.0*v[1] - 2.0*v[2]
824  -
825  4.0*e[0] - 4.0*e[1] - 4.0*e[2]
826  -
827  boxVol
828  );
829 
830  if (pt[0] <= centrePt[0])
831  {
832  if (pt[1] <= centrePt[1])
833  {
834  if (pt[2] <= centrePt[2])
835  {
836  vol = boxVol + v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
837  }
838  else
839  {
840  vol = v[2] + v[3] + e[1] + e[2];
841  }
842  }
843  else
844  {
845  if (pt[2] <= centrePt[2])
846  {
847  vol = v[1] + v[3] + e[0] + e[1];
848  }
849  else
850  {
851  vol = v[3] + e[1];
852  }
853  }
854  }
855  else
856  {
857  if (pt[1] <= centrePt[1])
858  {
859  if (pt[2] <= centrePt[2])
860  {
861  vol = v[0] + v[3] + e[0] + e[2];
862  }
863  else
864  {
865  vol = v[3] + e[2];
866  }
867  }
868  else
869  {
870  if (pt[2] <= centrePt[2])
871  {
872  vol = v[3] + e[0];
873  }
874  else
875  {
876  vol = v[3];
877  }
878  }
879  }
880  }
881 
882  return vol;
883  }
884 
885  template <int tmplDim, typename TmplVec>
886  double IntersectionVolCalculator<tmplDim,TmplVec>::getTwoPlaneVolume(
887  const Vec &pt,
888  const int orientDim
889  ) const
890  {
891  const int ZERO = orientDim;
892  const int ONE = (orientDim+1) % tmplDim;
893  const int TWO = (orientDim+2) % tmplDim;
894 
895  double vol = 0.0;
896  const double sphVol = getSphere().getVolume();
897  const Vec &centrePt = getSphere().getCentre();
898 
899  if ((square(pt[ONE]-centrePt[ONE]) + square(pt[TWO]-centrePt[TWO])) < square(getSphere().getRadius()))
900  {
901  Plane plane[tmplDim];
902  plane[ZERO] = Plane();
903  plane[ONE] = Plane(getNormal(ONE), pt);
904  plane[TWO] = Plane(getNormal(TWO), pt);
905 
906  const double halfSphVol = sphVol*0.5;
907  double s[tmplDim];
908  s[ZERO] = 0;
909  s[ONE] = getSphere().getSegmentVolume(plane[ONE]);
910  s[TWO] = getSphere().getSegmentVolume(plane[TWO]);
911  s[ONE] = ((s[ONE] > halfSphVol) ? (sphVol - s[ONE]) : s[ONE]);
912  s[TWO] = ((s[TWO] > halfSphVol) ? (sphVol - s[TWO]) : s[TWO]);
913 
914  Vec distVec(4.0*getSphere().getRadius());
915  distVec[ONE] = plane[ONE].getDistanceTo(centrePt);
916  distVec[TWO] = plane[TWO].getDistanceTo(centrePt);
917  const double coreVol =
918  2.0*getSphere().getVolume(
919  centrePt - Vec(distVec[ONE], distVec[TWO], distVec[ZERO]),
920  centrePt + Vec(distVec[ONE], distVec[TWO], distVec[ZERO])
921  );
922  double v[tmplDim];
923  v[ONE] = 0.50*(sphVol - 2.0*s[TWO] - coreVol);
924  v[TWO] = 0.50*(sphVol - 2.0*s[ONE] - coreVol);
925  v[ZERO] = 0.25*(sphVol - 2.0*v[ONE] - 2.0*v[TWO] - coreVol);
926 
927  if (pt[ONE] <= centrePt[ONE])
928  {
929  if (pt[TWO] <= centrePt[TWO])
930  {
931  vol = coreVol + v[ONE] + v[TWO] + v[ZERO];
932  }
933  else
934  {
935  vol = v[TWO] + v[ZERO];
936  }
937  }
938  else
939  {
940  if (pt[TWO] <= centrePt[TWO])
941  {
942  vol = v[ONE] + v[ZERO];
943  }
944  else
945  {
946  vol = v[ZERO];
947  }
948  }
949  }
950  else
951  {
952  if (pt[ONE] <= centrePt[ONE])
953  {
954  if (pt[TWO] <= centrePt[TWO])
955  {
956  vol =
957  sphVol
958  -
959  getSphere().getSegmentVolume(Plane(getNegNormal(ONE), pt))
960  -
961  getSphere().getSegmentVolume(Plane(getNegNormal(TWO), pt));
962  }
963  else
964  {
965  vol = getSphere().getSegmentVolume(Plane(getNormal(TWO), pt));
966  }
967  }
968  else
969  {
970  if (pt[TWO] <= centrePt[TWO])
971  {
972  vol = getSphere().getSegmentVolume(Plane(getNormal(ONE), pt));
973  }
974  else
975  {
976  vol = 0.0;
977  }
978  }
979  }
980 
981  return vol;
982  }
983 
984  template <int tmplDim, typename TmplVec>
985  double IntersectionVolCalculator<tmplDim,TmplVec>::getOutsidePointVolume(
986  const Vec &pt
987  ) const
988  {
989  double vol = 0.0;
990  const double sphVol = getSphere().getVolume();
991  const Vec &centrePt = getSphere().getCentre();
992  if (tmplDim == 2)
993  {
994  if (pt[0] <= centrePt[0])
995  {
996  if (pt[1] <= centrePt[1])
997  {
998  vol =
999  sphVol
1000  -
1001  getSphere().getSegmentVolume(Plane(getNegNormal(0), pt))
1002  -
1003  getSphere().getSegmentVolume(Plane(getNegNormal(1), pt));
1004  }
1005  else
1006  {
1007  vol = getSphere().getSegmentVolume(Plane(getNormal(1), pt));
1008  }
1009  }
1010  else
1011  {
1012  if (pt[1] <= centrePt[1])
1013  {
1014  vol = getSphere().getSegmentVolume(Plane(getNormal(0), pt));
1015  }
1016  else
1017  {
1018  vol = 0.0;
1019  }
1020  }
1021  }
1022  else if (tmplDim == 3)
1023  {
1024  const Vec diag = (centrePt - pt)*2.0;
1025  const Vec oppCorner = pt + diag;
1026  BasicBox box =
1027  BasicBox(
1028  componentMin(pt, oppCorner),
1029  componentMax(pt, oppCorner)
1030  );
1031 
1032  double s[tmplDim];
1033  double e[tmplDim];
1034  for (int i = 0; i < tmplDim; i++)
1035  {
1036  s[i] = getSphere().getSegmentVolume(Plane(getNormal(i), box.getMaxPt()));
1037  e[i] = getTwoPlaneVolume(box.getMaxPt(), i);
1038  }
1039  double v[tmplDim+1];
1040  v[0] = s[0] - 2.0*e[1] - 2.0*e[2];
1041  v[1] = s[1] - 2.0*e[0] - 2.0*e[2];
1042  v[2] = s[2] - 2.0*e[0] - 2.0*e[1];
1043  v[3] = sphVol - (4.0*e[0] + 4.0*e[1] + 4.0*e[2] + 2.0*v[0] + 2.0*v[1] + 2.0*v[2]);
1044 
1045  if (pt[0] <= centrePt[0])
1046  {
1047  if (pt[1] <= centrePt[1])
1048  {
1049  if (pt[2] <= centrePt[2])
1050  {
1051  vol = v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
1052  }
1053  else
1054  {
1055  vol = v[2] + e[0] + e[1];
1056  }
1057  }
1058  else
1059  {
1060  if (pt[2] <= centrePt[2])
1061  {
1062  vol = v[1] + e[0] + e[2];
1063  }
1064  else
1065  {
1066  vol = e[0];
1067  }
1068  }
1069  }
1070  else
1071  {
1072  if (pt[1] <= centrePt[1])
1073  {
1074  if (pt[2] <= centrePt[2])
1075  {
1076  vol = v[0] + e[1] + e[2];
1077  }
1078  else
1079  {
1080  vol = e[1];
1081  }
1082  }
1083  else
1084  {
1085  if (pt[2] <= centrePt[2])
1086  {
1087  vol = e[2];
1088  }
1089  else
1090  {
1091  vol = 0.0;
1092  }
1093  }
1094  }
1095  }
1096  return vol;
1097  }
1098 
1099 
1100  template <int tmplDim, typename TmplVec>
1101  double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
1102  const Vertex &vtx
1103  )
1104  {
1105  double vol = 0;
1106  if (getSphere().intersectsWith(vtx.getPoint()))
1107  {
1108  vol = getInsidePointVolume(vtx.getPoint());
1109  }
1110  else
1111  {
1112  vol = getOutsidePointVolume(vtx.getPoint());
1113  }
1114  return vol;
1115  }
1116 
1117  template <int tmplDim, typename TmplVec>
1118  double IntersectionVolCalculator<tmplDim,TmplVec>::getVertexVolume(
1119  const BasicSphere &sphere
1120  )
1121  {
1122  m_sphere = VolumeSphere(sphere);
1123  double vol = 0.0;
1124 
1125  double vtxVol[getVertexBox().getNumVertices()];
1126  for (int i = 0; i < getVertexBox().getNumVertices(); i++)
1127  {
1128  vtxVol[i] = getVolume(getVertexBox().getVertex(i));
1129  }
1130 
1131  if (tmplDim == 2)
1132  {
1133  vol = vtxVol[0] - vtxVol[1] - vtxVol[2] + vtxVol[3];
1134  }
1135  else if (tmplDim == 3)
1136  {
1137  vol =
1138  vtxVol[7] + vtxVol[6] + vtxVol[5] - vtxVol[4]
1139  -
1140  vtxVol[3] - vtxVol[2] - vtxVol[1] + vtxVol[0];
1141  }
1142 
1143  return vol;
1144  }
1145 
1146  template <int tmplDim, typename TmplVec>
1147  bool IntersectionVolCalculator<tmplDim,TmplVec>::sphereContainsBox(
1148  const BasicSphere &sphere
1149  ) const
1150  {
1151  for (int i = 0; i < getVertexBox().getNumVertices(); i++)
1152  {
1153  if (!sphere.intersectsWith(getVertexBox().getVertex(i).getPoint()))
1154  {
1155  return false;
1156  }
1157  }
1158  return true;
1159  }
1160 
1161  template <int tmplDim, typename TmplVec>
1162  double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
1163  const BasicSphere &sphere
1164  )
1165  {
1166  double vol = 0.0;
1167  if (getBox().intersectsWith(sphere))
1168  {
1169  if (sphereContainsBox(sphere))
1170  {
1171  vol = getBox().getVolume();
1172  }
1173  else if (getBox().contains(sphere))
1174  {
1175  vol = sphere.getVolume();
1176  }
1177  else
1178  {
1179  vol = getVertexVolume(sphere);
1180  }
1181  }
1182  return vol;
1183  }
1184  }
1185  }
1186 }