ESyS-Particle
4.0.1
|
00001 00002 // // 00003 // Copyright (c) 2003-2011 by The University of Queensland // 00004 // Earth Systems Science Computational Centre (ESSCC) // 00005 // http://www.uq.edu.au/esscc // 00006 // // 00007 // Primary Business: Brisbane, Queensland, Australia // 00008 // Licensed under the Open Software License version 3.0 // 00009 // http://www.opensource.org/licenses/osl-3.0.php // 00010 // // 00012 00013 00014 00015 #include <stdexcept> 00016 #include <sstream> 00017 00018 namespace esys 00019 { 00020 namespace lsm 00021 { 00022 namespace impl 00023 { 00024 inline double square(double val) 00025 { 00026 return val*val; 00027 } 00028 00029 template <int tmplDim, typename TmplVec> 00030 DimBasicBox<tmplDim,TmplVec>::DimBasicBox(const Vec &minPt, const Vec &maxPt) 00031 : m_minPt(minPt), 00032 m_maxPt(maxPt) 00033 { 00034 } 00035 00036 template <int tmplDim, typename TmplVec> 00037 const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMinPt() const 00038 { 00039 return m_minPt; 00040 } 00041 00042 template <int tmplDim, typename TmplVec> 00043 const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMaxPt() const 00044 { 00045 return m_maxPt; 00046 } 00047 00048 template <int tmplDim, typename TmplVec> 00049 double DimBasicBox<tmplDim,TmplVec>::getVolume() const 00050 { 00051 double product = 1.0; 00052 for (int i = 0; i < tmplDim; i++) 00053 { 00054 product *= (this->getMaxPt()[i] - this->getMinPt()[i]); 00055 } 00056 return product; 00057 } 00058 00059 template <int tmplDim, typename TmplVec> 00060 template <typename TmplSphere> 00061 bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const TmplSphere &sphere) const 00062 { 00063 double distSqrd = 0.0; 00064 for (int i = 0; i < tmplDim; i++) 00065 { 00066 if (sphere.getCentre()[i] < this->getMinPt()[i]) 00067 { 00068 distSqrd += square(sphere.getCentre()[i] - this->getMinPt()[i]); 00069 } 00070 else if (sphere.getCentre()[i] > this->getMaxPt()[i]) 00071 { 00072 distSqrd += square(sphere.getCentre()[i] - this->getMaxPt()[i]); 00073 } 00074 } 00075 return (distSqrd <= square(sphere.getRadius())); 00076 } 00077 00078 template <int tmplDim, typename TmplVec> 00079 bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const 00080 { 00081 for (int i = 0; (i < tmplDim); i++) 00082 { 00083 if ((this->getMinPt()[i] > pt[i]) || (pt[i] > this->getMaxPt()[i])) 00084 { 00085 return false; 00086 } 00087 } 00088 return true; 00089 } 00090 00091 template <int tmplDim, typename TmplVec> 00092 template <typename TmplSphere> 00093 bool DimBasicBox<tmplDim,TmplVec>::contains(const TmplSphere &sphere) const 00094 { 00095 for (int i = 0; (i < tmplDim); i++) 00096 { 00097 Vec pt = Vec(0.0); 00098 pt[i] = sphere.getRadius(); 00099 if 00100 ( 00101 !(intersectsWith(sphere.getCentre() + pt)) 00102 || 00103 !(intersectsWith(sphere.getCentre() - pt)) 00104 ) 00105 { 00106 return false; 00107 } 00108 } 00109 return true; 00110 } 00111 00112 template <int tmplDim, typename TmplVec> 00113 double DimPlane<tmplDim,TmplVec>::norm(const Vec &pt) 00114 { 00115 double sum = 0.0; 00116 for (int i = 0; i < tmplDim; i++) 00117 { 00118 sum += pt[i]*pt[i]; 00119 } 00120 return sqrt(sum); 00121 } 00122 00123 template <int tmplDim, typename TmplVec> 00124 double DimPlane<tmplDim,TmplVec>::dot(const Vec &p1, const Vec &p2) 00125 { 00126 double sum = 0.0; 00127 for (int i = 0; i < tmplDim; i++) 00128 { 00129 sum += p1[i]*p2[i]; 00130 } 00131 return sum; 00132 } 00133 00134 template <int tmplDim, typename TmplVec> 00135 DimPlane<tmplDim,TmplVec>::DimPlane() : m_normal(), m_pt(), m_invNormalNorm(0.0) 00136 { 00137 } 00138 00139 template <int tmplDim, typename TmplVec> 00140 DimPlane<tmplDim,TmplVec>::DimPlane(const Vec &normal, const Vec &pt) 00141 : m_normal(normal), 00142 m_pt(pt), 00143 m_invNormalNorm((1.0/norm(normal))) 00144 { 00145 } 00146 00147 template <int tmplDim, typename TmplVec> 00148 DimPlane<tmplDim,TmplVec>::DimPlane(const DimPlane &plane) 00149 : m_normal(plane.m_normal), 00150 m_pt(plane.m_pt), 00151 m_invNormalNorm(plane.m_invNormalNorm) 00152 { 00153 } 00154 00155 template <int tmplDim, typename TmplVec> 00156 DimPlane<tmplDim,TmplVec> &DimPlane<tmplDim,TmplVec>::operator=(const DimPlane &plane) 00157 { 00158 m_normal = plane.m_normal; 00159 m_pt = plane.m_pt; 00160 m_invNormalNorm = plane.m_invNormalNorm; 00161 return *this; 00162 } 00163 00164 template <int tmplDim, typename TmplVec> 00165 double DimPlane<tmplDim,TmplVec>::getSignedDistanceTo(const Vec &pt) const 00166 { 00167 // http://mathworld.wolfram.com/Point-PlaneDistance.html 00168 return 00169 ( 00170 (dot(m_normal, pt) - dot(m_normal, m_pt))*m_invNormalNorm 00171 ); 00172 } 00173 00174 template <int tmplDim, typename TmplVec> 00175 double DimPlane<tmplDim,TmplVec>::getDistanceTo(const Vec &pt) const 00176 { 00177 return fabs(getSignedDistanceTo(pt)); 00178 } 00179 00180 template <int tmplDim, typename TmplVec> 00181 const typename DimPlane<tmplDim,TmplVec>::Vec &DimPlane<tmplDim,TmplVec>::getNormal() const 00182 { 00183 return m_normal; 00184 } 00185 00186 template <int tmplDim, typename TmplVec> 00187 const double DimBasicSphere<tmplDim,TmplVec>::FOUR_THIRDS_PI = (4.0/3.0)*M_PI; 00188 00189 template <int tmplDim, typename TmplVec> 00190 const double DimBasicSphere<tmplDim,TmplVec>::ONE_THIRD_PI = M_PI/3.0; 00191 00192 template <int tmplDim, typename TmplVec> 00193 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere() 00194 : m_centre(), 00195 m_radius(0.0) 00196 { 00197 } 00198 00199 template <int tmplDim, typename TmplVec> 00200 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const Vec ¢rePt, double radius) 00201 : m_centre(centrePt), 00202 m_radius(radius) 00203 { 00204 } 00205 00206 template <int tmplDim, typename TmplVec> 00207 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const DimBasicSphere &sphere) 00208 : m_centre(sphere.m_centre), 00209 m_radius(sphere.m_radius) 00210 { 00211 } 00212 00213 template <int tmplDim, typename TmplVec> 00214 DimBasicSphere<tmplDim,TmplVec> &DimBasicSphere<tmplDim,TmplVec>::operator=(const DimBasicSphere &sphere) 00215 { 00216 m_centre = sphere.m_centre; 00217 m_radius = sphere.m_radius; 00218 return *this; 00219 } 00220 00221 template <int tmplDim, typename TmplVec> 00222 double DimBasicSphere<tmplDim,TmplVec>::getRadius() const 00223 { 00224 return m_radius; 00225 } 00226 00227 template <int tmplDim, typename TmplVec> 00228 const typename DimBasicSphere<tmplDim,TmplVec>::Vec & 00229 DimBasicSphere<tmplDim,TmplVec>::getCentre() const 00230 { 00231 return m_centre; 00232 } 00233 00234 template <int tmplDim, typename TmplVec> 00235 double DimBasicSphere<tmplDim,TmplVec>::getVolume() const 00236 { 00237 return (tmplDim == 2) ? M_PI*getRadius()*getRadius() : FOUR_THIRDS_PI*getRadius()*getRadius()*getRadius(); 00238 } 00239 00240 inline void checkDomain(double r, double x1, double y1, double x2, double y2) 00241 { 00242 const double rSqrd = r*r; 00243 const double x1Sqrd = x1*x1; 00244 const double x2Sqrd = x2*x2; 00245 const double y1Sqrd = y1*y1; 00246 const double y2Sqrd = y2*y2; 00247 if 00248 ( 00249 ((rSqrd - x1Sqrd - y1Sqrd) < 0) 00250 || 00251 ((rSqrd - x1Sqrd - y2Sqrd) < 0) 00252 || 00253 ((rSqrd - x2Sqrd - y1Sqrd) < 0) 00254 || 00255 ((rSqrd - x2Sqrd - y2Sqrd) < 0) 00256 ) 00257 { 00258 std::stringstream msg; 00259 msg 00260 << "Invalid rectangular domain for sphere integration, points in domain " 00261 << "(" << x1 << "," << y1 << "), (" << x2 << "," << y2 << " lie outside " 00262 << "sphere of radius " << r << " centred at the origin."; 00263 throw std::runtime_error(msg.str()); 00264 } 00265 } 00266 00267 template <int tmplDim, typename TmplVec> 00268 double DimBasicSphere<tmplDim,TmplVec>::getVolume(const Vec &minPt, const Vec &maxPt, const int dimX, const int dimY) const 00269 { 00270 double vol = 0.0; 00271 if ((tmplDim == 2) || (tmplDim == 3)) 00272 { 00273 if (minPt[dimX] != maxPt[dimX]) 00274 { 00275 const double x1 = minPt[dimX] - getCentre()[dimX]; 00276 const double x2 = maxPt[dimX] - getCentre()[dimX]; 00277 const double r = getRadius(); 00278 00279 if (tmplDim == 2) 00280 { 00281 const double rSqrd = r*r; 00282 const double x1Sqrd = x1*x1; 00283 const double x2Sqrd = x2*x2; 00284 00285 vol = 00286 ( 00287 rSqrd*asin(x2/r) 00288 + 00289 x2*sqrt(rSqrd-x2Sqrd) 00290 - 00291 rSqrd*asin(x1/r) 00292 - 00293 x1*sqrt(rSqrd-x1Sqrd) 00294 )*0.5; 00295 } 00296 else if (tmplDim == 3) 00297 { 00298 if (minPt[dimY] != maxPt[dimY]) 00299 { 00300 const double y1 = minPt[dimY] - getCentre()[dimY]; 00301 const double y2 = maxPt[dimY] - getCentre()[dimY]; 00302 00303 checkDomain(r, x1, y1, x2, y2); 00304 00305 // 00306 // Matlab6/maple generated code: 00307 // 00308 // syms x y x1 x2 y1 y2 r real 00309 // sphereIntegral = int(int('sqrt(r^2-x^2-y^2)', x, x1, x2), y, y1, y2) 00310 // maple('readlib(codegen)') 00311 // maple('readlib(optimize)') 00312 // optimizedIntegral = maple('optimize', sphereIntegral, 'tryhard') 00313 // maple('cost', optimizedIntegral) 00314 // 00315 //43*additions+82*multiplications+6*divisions+22*functions+42*assignments 00316 00317 const double t30 = y2*y2; //y2^2; 00318 const double t31 = x2*x2; //x2^2; 00319 const double t36 = r*r; //r^2; 00320 const double t59 = t31-t36; 00321 const double t40 = sqrt(-t30-t59); //pow(-t30-t59,1/2); 00322 const double t10 = 1.0/t40; 00323 const double t32 = x1*x1; //x1^2; 00324 const double t54 = t32-t36; 00325 const double t42 = sqrt(-t30-t54); //pow(-t30-t54,1/2); 00326 const double t14 = 1.0/t42; 00327 const double t64 = -atan(t10*x2)+atan(t14*x1); 00328 const double t27 = y1*y1; //y1^2; 00329 const double t39 = sqrt(-t27-t59); //pow(-t27-t59,1/2); 00330 const double t9 = 1.0/t39; 00331 const double t41 = sqrt(-t27-t54); //pow(-t27-t54,1/2); 00332 const double t12 = 1.0/t41; 00333 const double t63 = atan(t12*x1)-atan(t9*x2); 00334 const double t62 = -atan(y2*t14)+atan(y1*t12); 00335 const double t61 = -atan(y1*t9)+atan(t10*y2); 00336 const double t37 = sqrt(t31); //pow(t31,1/2); 00337 const double t21 = 1.0/t37; 00338 const double t60 = t21*t9; 00339 const double t38 = sqrt(t32); //pow(t32,1/2); 00340 const double t24 = 1.0/t38; 00341 const double t58 = t24*t14; 00342 const double t57 = t24*t12; 00343 const double t56 = t37*t38; 00344 const double t55 = t21*t10; 00345 const double t53 = 2.0*x2; 00346 const double t52 = 2.0*x1; 00347 const double t51 = t42*t56; 00348 const double t28 = t27*y1; 00349 const double t50 = t28-t36*y1; 00350 const double t34 = t30*y2; 00351 const double t49 = t34-t36*y2; 00352 const double t48 = t41*t51; 00353 const double t35 = t36*r; 00354 const double t33 = t31*x2; 00355 const double t29 = t32*x1; 00356 const double t26 = r*y2; 00357 const double t25 = y1*r; 00358 vol = 00359 (-1.0/6.0) 00360 * 00361 ( 00362 (-2.0*t33*y1-t50*t53)*t40*t48 00363 + 00364 ( 00365 (2.0*t33*y2+t49*t53)*t48 00366 + 00367 ( 00368 (2.0*t29*y1+t50*t52)*t51 00369 + 00370 ( 00371 (-2.0*t29*y2-t49*t52)*t56 00372 + 00373 ( 00374 ( 00375 atan((t25+t54)*t57) 00376 + 00377 atan((-t26+t54)*t58) 00378 - 00379 atan((t26+t54)*t58) 00380 - 00381 atan((-t25+t54)*t57) 00382 )*t35*t37*x1 00383 + 00384 ( 00385 ( 00386 -atan((-t26+t59)*t55) 00387 - 00388 atan((t25+t59)*t60) 00389 + 00390 atan((-t25+t59)*t60) 00391 + 00392 atan((t26+t59)*t55) 00393 )*t35*x2 00394 + 00395 (-t64*t34+t61*t33+t62*t29+t63*t28+3.0*(t64*y2-t63*y1-t61*x2-t62*x1)*t36)*t37 00396 )*t38 00397 )*t42 00398 )*t41 00399 )*t40 00400 )*t39 00401 )*t14*t9*t55*t57; 00402 } 00403 } 00404 } 00405 } 00406 00407 return vol; 00408 } 00409 00410 template <int tmplDim, typename TmplVec> 00411 bool DimBasicSphere<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const 00412 { 00413 double distSqrd = 0.0; 00414 for (int i = 0; i < tmplDim; i++) 00415 { 00416 distSqrd += square(getCentre()[i] - pt[i]); 00417 } 00418 return (distSqrd <= square(getRadius())); 00419 } 00420 00421 template <int tmplDim, typename TmplVec> 00422 double DimBasicSphere<tmplDim,TmplVec>::getSegmentVolume(const Plane &plane) const 00423 { 00424 double vol = 0.0; 00425 if ((tmplDim == 2) || (tmplDim == 3)) 00426 { 00427 const double signedD = plane.getSignedDistanceTo(getCentre()); 00428 const double d = fabs(signedD); 00429 if (d < getRadius()) 00430 { 00431 if (tmplDim == 2) 00432 { 00433 // http://mathworld.wolfram.com/CircularSegment.html 00434 const double rSqrd = getRadius()*getRadius(); 00435 vol = rSqrd*acos(d/getRadius()) - d*sqrt(rSqrd - d*d); 00436 } 00437 else if (tmplDim == 3) 00438 { 00439 // http://mathworld.wolfram.com/SphericalCap.html 00440 const double h = getRadius() - d; 00441 vol = ONE_THIRD_PI*h*h*(3.0*getRadius()-h); 00442 } 00443 vol = ((signedD < 0) ? vol : getVolume() - vol); 00444 } 00445 } 00446 return vol; 00447 } 00448 00449 template <int tmplDim, typename TmplVec> 00450 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec 00451 IntersectionVolCalculator<tmplDim,TmplVec>::getNormal(int dim) 00452 { 00453 Vec n = Vec(0.0); 00454 n[dim] = 1.0; 00455 return n; 00456 } 00457 00458 template <int tmplDim, typename TmplVec> 00459 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec 00460 IntersectionVolCalculator<tmplDim,TmplVec>::getNegNormal(int dim) 00461 { 00462 Vec n = Vec(0.0); 00463 n[dim] = -1.0; 00464 return n; 00465 } 00466 00467 template <int tmplDim, typename TmplVec> 00468 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere() 00469 : m_sphere(), 00470 m_volume(0.0) 00471 { 00472 } 00473 00474 template <int tmplDim, typename TmplVec> 00475 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere( 00476 const BasicSphere &sphere 00477 ) 00478 : m_sphere(sphere), 00479 m_volume(sphere.getVolume()) 00480 { 00481 } 00482 00483 template <int tmplDim, typename TmplVec> 00484 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere( 00485 const VolumeSphere &sphere 00486 ) 00487 : m_sphere(BasicSphere(sphere.getCentre(), sphere.getRadius())), 00488 m_volume(sphere.m_volume) 00489 { 00490 } 00491 00492 template <int tmplDim, typename TmplVec> 00493 typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere & 00494 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::operator=( 00495 const VolumeSphere &sphere 00496 ) 00497 { 00498 m_sphere = sphere.m_sphere; 00499 m_volume = sphere.m_volume; 00500 return *this; 00501 } 00502 00503 template <int tmplDim, typename TmplVec> 00504 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getRadius() const 00505 { 00506 return m_sphere.getRadius(); 00507 } 00508 00509 template <int tmplDim, typename TmplVec> 00510 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec & 00511 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getCentre() const 00512 { 00513 return m_sphere.getCentre(); 00514 } 00515 00516 template <int tmplDim, typename TmplVec> 00517 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume() const 00518 { 00519 return m_volume; 00520 } 00521 00522 template <int tmplDim, typename TmplVec> 00523 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume( 00524 const Vec &minPt, 00525 const Vec &maxPt, 00526 const int dimX, 00527 const int dimY 00528 ) const 00529 { 00530 return m_sphere.getVolume(minPt, maxPt, dimX, dimY); 00531 } 00532 00533 template <int tmplDim, typename TmplVec> 00534 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::calcVolume() const 00535 { 00536 return m_sphere.getVolume(); 00537 } 00538 00539 template <int tmplDim, typename TmplVec> 00540 bool IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::intersectsWith(const Vec &pt) const 00541 { 00542 return m_sphere.intersectsWith(pt); 00543 } 00544 00545 template <int tmplDim, typename TmplVec> 00546 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getSegmentVolume(const Plane &plane) const 00547 { 00548 return m_sphere.getSegmentVolume(plane); 00549 } 00550 00551 template <int tmplDim, typename TmplVec> 00552 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex() : m_pt() 00553 { 00554 } 00555 00556 template <int tmplDim, typename TmplVec> 00557 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vec &pt) : m_pt(pt) 00558 { 00559 } 00560 00561 template <int tmplDim, typename TmplVec> 00562 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vertex &vtx) : m_pt(vtx.m_pt) 00563 { 00564 } 00565 00566 template <int tmplDim, typename TmplVec> 00567 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex & 00568 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::operator=(const Vertex &vtx) 00569 { 00570 m_pt = vtx.m_pt; 00571 return *this; 00572 } 00573 00574 template <int tmplDim, typename TmplVec> 00575 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec & 00576 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::getPoint() const 00577 { 00578 return m_pt; 00579 } 00580 00581 template <int tmplDim, typename TmplVec> 00582 void IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::setPoint(const Vec &pt) 00583 { 00584 m_pt = pt; 00585 } 00586 00587 template <int tmplDim, typename TmplVec> 00588 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox( 00589 const BasicBox &box 00590 ) 00591 : BasicBox(box), 00592 m_vertexArray() 00593 { 00594 createVertices(); 00595 } 00596 00597 template <int tmplDim, typename TmplVec> 00598 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox( 00599 const VertexBox &box 00600 ) 00601 : BasicBox(box) 00602 { 00603 for (int i = 0; i < getNumVertices(); i++) 00604 { 00605 m_vertexArray[i] = box.m_vertexArray[i]; 00606 } 00607 } 00608 00609 template <int tmplDim, typename TmplVec> 00610 typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox & 00611 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::operator=( 00612 const VertexBox &box 00613 ) 00614 { 00615 BasicBox::operator=(box); 00616 for (int i = 0; i < getNumVertices(); i++) 00617 { 00618 m_vertexArray[i] = box.m_vertexArray[i]; 00619 } 00620 return *this; 00621 } 00622 00623 template <int tmplDim, typename TmplVec> 00624 void IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::createVertices() 00625 { 00626 int j = 0; 00627 m_vertexArray[j].setPoint(this->getMinPt()); 00628 int i = 0; 00629 for (j++; i < tmplDim; i++, j++) 00630 { 00631 Vec pt = this->getMinPt(); 00632 pt[i] = this->getMaxPt()[i]; 00633 m_vertexArray[j].setPoint(pt); 00634 } 00635 00636 m_vertexArray[j].setPoint(this->getMaxPt()); 00637 for (i = 0, j++; i < tmplDim && j < s_numVertices; i++, j++) 00638 { 00639 Vec pt = this->getMaxPt(); 00640 pt[i] = this->getMinPt()[i]; 00641 m_vertexArray[j] = pt; 00642 } 00643 } 00644 00645 template <int tmplDim, typename TmplVec> 00646 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex & 00647 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getVertex(int i) const 00648 { 00649 return m_vertexArray[i]; 00650 } 00651 00652 template <int tmplDim, typename TmplVec> 00653 int IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getNumVertices() 00654 { 00655 return s_numVertices; 00656 } 00657 00658 template <int tmplDim, typename TmplVec> 00659 IntersectionVolCalculator<tmplDim,TmplVec>::IntersectionVolCalculator( 00660 const BasicBox &box 00661 ) 00662 : m_sphere(BasicSphere(Vec(), 1.0)), 00663 m_box(box) 00664 { 00665 } 00666 00667 template <int tmplDim, typename TmplVec> 00668 const typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere & 00669 IntersectionVolCalculator<tmplDim,TmplVec>::getSphere() const 00670 { 00671 return m_sphere; 00672 } 00673 00674 template <int tmplDim, typename TmplVec> 00675 void IntersectionVolCalculator<tmplDim,TmplVec>::setSphere( 00676 const BasicSphere &sphere 00677 ) 00678 { 00679 m_sphere = sphere; 00680 } 00681 00682 template <int tmplDim, typename TmplVec> 00683 const typename IntersectionVolCalculator<tmplDim,TmplVec>::BasicBox & 00684 IntersectionVolCalculator<tmplDim,TmplVec>::getBox() const 00685 { 00686 return m_box; 00687 } 00688 00689 template <int tmplDim, typename TmplVec> 00690 const typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox & 00691 IntersectionVolCalculator<tmplDim,TmplVec>::getVertexBox() const 00692 { 00693 return m_box; 00694 } 00695 00696 template <int tmplDim, typename TmplVec> 00697 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec 00698 IntersectionVolCalculator<tmplDim,TmplVec>::componentMin( 00699 const Vec &p1, 00700 const Vec &p2 00701 ) 00702 { 00703 Vec m; 00704 for (int i = 0; i < tmplDim; i++) 00705 { 00706 m[i] = min(p1[i], p2[i]); 00707 } 00708 return m; 00709 } 00710 00711 template <int tmplDim, typename TmplVec> 00712 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec 00713 IntersectionVolCalculator<tmplDim,TmplVec>::componentMax( 00714 const Vec &p1, 00715 const Vec &p2 00716 ) 00717 { 00718 Vec m; 00719 for (int i = 0; i < tmplDim; i++) 00720 { 00721 m[i] = max(p1[i], p2[i]); 00722 } 00723 return m; 00724 } 00725 00726 template <int tmplDim, typename TmplVec> 00727 double IntersectionVolCalculator<tmplDim,TmplVec>::getInsidePointVolume( 00728 const Vec &pt 00729 ) const 00730 { 00731 double vol = 0.0; 00732 const Vec ¢rePt = getSphere().getCentre(); 00733 const Vec diag = (getSphere().getCentre() - pt)*2.0; 00734 const Vec oppCorner = pt + diag; 00735 BasicBox box = 00736 BasicBox( 00737 componentMin(pt, oppCorner), 00738 componentMax(pt, oppCorner) 00739 ); 00740 const double boxVol = box.getVolume(); 00741 const double sphVol = getSphere().getVolume(); 00742 00743 double s[tmplDim]; 00744 double v[tmplDim+1]; 00745 for (int i = 0; i < tmplDim; i++) 00746 { 00747 s[i] = getSphere().getSegmentVolume(Plane(getNormal((i+1)%tmplDim), box.getMaxPt())); 00748 } 00749 if (tmplDim == 2) 00750 { 00751 v[0] = 0.50*(sphVol - 2.0*s[0] - boxVol); 00752 v[1] = 0.50*(sphVol - 2.0*s[1] - boxVol); 00753 v[2] = 0.25*(sphVol - 2.0*v[0] - 2.0*v[1] - boxVol); 00754 00755 if (pt[0] <= centrePt[0]) 00756 { 00757 if (pt[1] <= centrePt[1]) 00758 { 00759 vol = boxVol + v[0] + v[1] + v[2]; 00760 } 00761 else 00762 { 00763 vol = v[1] + v[2]; 00764 } 00765 } 00766 else 00767 { 00768 if (pt[1] <= centrePt[1]) 00769 { 00770 vol = v[0] + v[2]; 00771 } 00772 else 00773 { 00774 vol = v[2]; 00775 } 00776 } 00777 } 00778 else if (tmplDim == 3) 00779 { 00780 v[0] = 00781 0.500*( 00782 2.0*getSphere().getVolume( 00783 box.getMinPt(), 00784 box.getMaxPt(), 00785 1, 00786 2 00787 ) 00788 - 00789 boxVol 00790 ); 00791 v[1] = 00792 0.500*( 00793 2.0*getSphere().getVolume( 00794 box.getMinPt(), 00795 box.getMaxPt(), 00796 0, 00797 2 00798 ) 00799 - 00800 boxVol 00801 ); 00802 v[2] = 00803 0.500*( 00804 2.0*getSphere().getVolume( 00805 box.getMinPt(), 00806 box.getMaxPt(), 00807 0, 00808 1 00809 ) 00810 - 00811 boxVol 00812 ); 00813 00814 double e[3]; 00815 e[0] = 0.250*(sphVol - 2.0*s[1] - 2.0*v[0] - 2.0*v[1] - boxVol); 00816 e[1] = 0.250*(sphVol - 2.0*s[2] - 2.0*v[1] - 2.0*v[2] - boxVol); 00817 e[2] = 0.250*(sphVol - 2.0*s[0] - 2.0*v[0] - 2.0*v[2] - boxVol); 00818 00819 v[3] = 00820 0.125*( 00821 sphVol 00822 - 00823 2.0*v[0] - 2.0*v[1] - 2.0*v[2] 00824 - 00825 4.0*e[0] - 4.0*e[1] - 4.0*e[2] 00826 - 00827 boxVol 00828 ); 00829 00830 if (pt[0] <= centrePt[0]) 00831 { 00832 if (pt[1] <= centrePt[1]) 00833 { 00834 if (pt[2] <= centrePt[2]) 00835 { 00836 vol = boxVol + v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2]; 00837 } 00838 else 00839 { 00840 vol = v[2] + v[3] + e[1] + e[2]; 00841 } 00842 } 00843 else 00844 { 00845 if (pt[2] <= centrePt[2]) 00846 { 00847 vol = v[1] + v[3] + e[0] + e[1]; 00848 } 00849 else 00850 { 00851 vol = v[3] + e[1]; 00852 } 00853 } 00854 } 00855 else 00856 { 00857 if (pt[1] <= centrePt[1]) 00858 { 00859 if (pt[2] <= centrePt[2]) 00860 { 00861 vol = v[0] + v[3] + e[0] + e[2]; 00862 } 00863 else 00864 { 00865 vol = v[3] + e[2]; 00866 } 00867 } 00868 else 00869 { 00870 if (pt[2] <= centrePt[2]) 00871 { 00872 vol = v[3] + e[0]; 00873 } 00874 else 00875 { 00876 vol = v[3]; 00877 } 00878 } 00879 } 00880 } 00881 00882 return vol; 00883 } 00884 00885 template <int tmplDim, typename TmplVec> 00886 double IntersectionVolCalculator<tmplDim,TmplVec>::getTwoPlaneVolume( 00887 const Vec &pt, 00888 const int orientDim 00889 ) const 00890 { 00891 const int ZERO = orientDim; 00892 const int ONE = (orientDim+1) % tmplDim; 00893 const int TWO = (orientDim+2) % tmplDim; 00894 00895 double vol = 0.0; 00896 const double sphVol = getSphere().getVolume(); 00897 const Vec ¢rePt = getSphere().getCentre(); 00898 00899 if ((square(pt[ONE]-centrePt[ONE]) + square(pt[TWO]-centrePt[TWO])) < square(getSphere().getRadius())) 00900 { 00901 Plane plane[tmplDim]; 00902 plane[ZERO] = Plane(); 00903 plane[ONE] = Plane(getNormal(ONE), pt); 00904 plane[TWO] = Plane(getNormal(TWO), pt); 00905 00906 const double halfSphVol = sphVol*0.5; 00907 double s[tmplDim]; 00908 s[ZERO] = 0; 00909 s[ONE] = getSphere().getSegmentVolume(plane[ONE]); 00910 s[TWO] = getSphere().getSegmentVolume(plane[TWO]); 00911 s[ONE] = ((s[ONE] > halfSphVol) ? (sphVol - s[ONE]) : s[ONE]); 00912 s[TWO] = ((s[TWO] > halfSphVol) ? (sphVol - s[TWO]) : s[TWO]); 00913 00914 Vec distVec(4.0*getSphere().getRadius()); 00915 distVec[ONE] = plane[ONE].getDistanceTo(centrePt); 00916 distVec[TWO] = plane[TWO].getDistanceTo(centrePt); 00917 const double coreVol = 00918 2.0*getSphere().getVolume( 00919 centrePt - Vec(distVec[ONE], distVec[TWO], distVec[ZERO]), 00920 centrePt + Vec(distVec[ONE], distVec[TWO], distVec[ZERO]) 00921 ); 00922 double v[tmplDim]; 00923 v[ONE] = 0.50*(sphVol - 2.0*s[TWO] - coreVol); 00924 v[TWO] = 0.50*(sphVol - 2.0*s[ONE] - coreVol); 00925 v[ZERO] = 0.25*(sphVol - 2.0*v[ONE] - 2.0*v[TWO] - coreVol); 00926 00927 if (pt[ONE] <= centrePt[ONE]) 00928 { 00929 if (pt[TWO] <= centrePt[TWO]) 00930 { 00931 vol = coreVol + v[ONE] + v[TWO] + v[ZERO]; 00932 } 00933 else 00934 { 00935 vol = v[TWO] + v[ZERO]; 00936 } 00937 } 00938 else 00939 { 00940 if (pt[TWO] <= centrePt[TWO]) 00941 { 00942 vol = v[ONE] + v[ZERO]; 00943 } 00944 else 00945 { 00946 vol = v[ZERO]; 00947 } 00948 } 00949 } 00950 else 00951 { 00952 if (pt[ONE] <= centrePt[ONE]) 00953 { 00954 if (pt[TWO] <= centrePt[TWO]) 00955 { 00956 vol = 00957 sphVol 00958 - 00959 getSphere().getSegmentVolume(Plane(getNegNormal(ONE), pt)) 00960 - 00961 getSphere().getSegmentVolume(Plane(getNegNormal(TWO), pt)); 00962 } 00963 else 00964 { 00965 vol = getSphere().getSegmentVolume(Plane(getNormal(TWO), pt)); 00966 } 00967 } 00968 else 00969 { 00970 if (pt[TWO] <= centrePt[TWO]) 00971 { 00972 vol = getSphere().getSegmentVolume(Plane(getNormal(ONE), pt)); 00973 } 00974 else 00975 { 00976 vol = 0.0; 00977 } 00978 } 00979 } 00980 00981 return vol; 00982 } 00983 00984 template <int tmplDim, typename TmplVec> 00985 double IntersectionVolCalculator<tmplDim,TmplVec>::getOutsidePointVolume( 00986 const Vec &pt 00987 ) const 00988 { 00989 double vol = 0.0; 00990 const double sphVol = getSphere().getVolume(); 00991 const Vec ¢rePt = getSphere().getCentre(); 00992 if (tmplDim == 2) 00993 { 00994 if (pt[0] <= centrePt[0]) 00995 { 00996 if (pt[1] <= centrePt[1]) 00997 { 00998 vol = 00999 sphVol 01000 - 01001 getSphere().getSegmentVolume(Plane(getNegNormal(0), pt)) 01002 - 01003 getSphere().getSegmentVolume(Plane(getNegNormal(1), pt)); 01004 } 01005 else 01006 { 01007 vol = getSphere().getSegmentVolume(Plane(getNormal(1), pt)); 01008 } 01009 } 01010 else 01011 { 01012 if (pt[1] <= centrePt[1]) 01013 { 01014 vol = getSphere().getSegmentVolume(Plane(getNormal(0), pt)); 01015 } 01016 else 01017 { 01018 vol = 0.0; 01019 } 01020 } 01021 } 01022 else if (tmplDim == 3) 01023 { 01024 const Vec diag = (centrePt - pt)*2.0; 01025 const Vec oppCorner = pt + diag; 01026 BasicBox box = 01027 BasicBox( 01028 componentMin(pt, oppCorner), 01029 componentMax(pt, oppCorner) 01030 ); 01031 01032 double s[tmplDim]; 01033 double e[tmplDim]; 01034 for (int i = 0; i < tmplDim; i++) 01035 { 01036 s[i] = getSphere().getSegmentVolume(Plane(getNormal(i), box.getMaxPt())); 01037 e[i] = getTwoPlaneVolume(box.getMaxPt(), i); 01038 } 01039 double v[tmplDim+1]; 01040 v[0] = s[0] - 2.0*e[1] - 2.0*e[2]; 01041 v[1] = s[1] - 2.0*e[0] - 2.0*e[2]; 01042 v[2] = s[2] - 2.0*e[0] - 2.0*e[1]; 01043 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]); 01044 01045 if (pt[0] <= centrePt[0]) 01046 { 01047 if (pt[1] <= centrePt[1]) 01048 { 01049 if (pt[2] <= centrePt[2]) 01050 { 01051 vol = v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2]; 01052 } 01053 else 01054 { 01055 vol = v[2] + e[0] + e[1]; 01056 } 01057 } 01058 else 01059 { 01060 if (pt[2] <= centrePt[2]) 01061 { 01062 vol = v[1] + e[0] + e[2]; 01063 } 01064 else 01065 { 01066 vol = e[0]; 01067 } 01068 } 01069 } 01070 else 01071 { 01072 if (pt[1] <= centrePt[1]) 01073 { 01074 if (pt[2] <= centrePt[2]) 01075 { 01076 vol = v[0] + e[1] + e[2]; 01077 } 01078 else 01079 { 01080 vol = e[1]; 01081 } 01082 } 01083 else 01084 { 01085 if (pt[2] <= centrePt[2]) 01086 { 01087 vol = e[2]; 01088 } 01089 else 01090 { 01091 vol = 0.0; 01092 } 01093 } 01094 } 01095 } 01096 return vol; 01097 } 01098 01099 01100 template <int tmplDim, typename TmplVec> 01101 double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume( 01102 const Vertex &vtx 01103 ) 01104 { 01105 double vol = 0; 01106 if (getSphere().intersectsWith(vtx.getPoint())) 01107 { 01108 vol = getInsidePointVolume(vtx.getPoint()); 01109 } 01110 else 01111 { 01112 vol = getOutsidePointVolume(vtx.getPoint()); 01113 } 01114 return vol; 01115 } 01116 01117 template <int tmplDim, typename TmplVec> 01118 double IntersectionVolCalculator<tmplDim,TmplVec>::getVertexVolume( 01119 const BasicSphere &sphere 01120 ) 01121 { 01122 m_sphere = VolumeSphere(sphere); 01123 double vol = 0.0; 01124 01125 double vtxVol[getVertexBox().getNumVertices()]; 01126 for (int i = 0; i < getVertexBox().getNumVertices(); i++) 01127 { 01128 vtxVol[i] = getVolume(getVertexBox().getVertex(i)); 01129 } 01130 01131 if (tmplDim == 2) 01132 { 01133 vol = vtxVol[0] - vtxVol[1] - vtxVol[2] + vtxVol[3]; 01134 } 01135 else if (tmplDim == 3) 01136 { 01137 vol = 01138 vtxVol[7] + vtxVol[6] + vtxVol[5] - vtxVol[4] 01139 - 01140 vtxVol[3] - vtxVol[2] - vtxVol[1] + vtxVol[0]; 01141 } 01142 01143 return vol; 01144 } 01145 01146 template <int tmplDim, typename TmplVec> 01147 bool IntersectionVolCalculator<tmplDim,TmplVec>::sphereContainsBox( 01148 const BasicSphere &sphere 01149 ) const 01150 { 01151 for (int i = 0; i < getVertexBox().getNumVertices(); i++) 01152 { 01153 if (!sphere.intersectsWith(getVertexBox().getVertex(i).getPoint())) 01154 { 01155 return false; 01156 } 01157 } 01158 return true; 01159 } 01160 01161 template <int tmplDim, typename TmplVec> 01162 double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume( 01163 const BasicSphere &sphere 01164 ) 01165 { 01166 double vol = 0.0; 01167 if (getBox().intersectsWith(sphere)) 01168 { 01169 if (sphereContainsBox(sphere)) 01170 { 01171 vol = getBox().getVolume(); 01172 } 01173 else if (getBox().contains(sphere)) 01174 { 01175 vol = sphere.getVolume(); 01176 } 01177 else 01178 { 01179 vol = getVertexVolume(sphere); 01180 } 01181 } 01182 return vol; 01183 } 01184 } 01185 } 01186 }