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