cgma
volumes.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <map>
00003 #include <cfloat>
00004 
00005 #include <cassert>
00006 
00007 #include "MCNPInput.hpp"
00008 #include "volumes.hpp"
00009 #include "geometry.hpp"
00010 #include "options.hpp"
00011 
00012 
00013 static Vector3d origin(0,0,0);
00014 
00015 
00016 iBase_EntityHandle makeWorldSphere( iGeom_Instance& igm, double world_size ){
00017   iBase_EntityHandle world_sphere;
00018   int igm_result;
00019   // Note: I tried using createBrick instead of createSphere to bound the universe with a box
00020   // instead of a sphere.  This worked but led to a substantial increase in run times and
00021   // memory usage, so should be avoided.
00022   iGeom_createSphere( igm, world_size, &world_sphere, &igm_result);
00023   CHECK_IGEOM( igm_result, "making world sphere" );
00024   return world_sphere;
00025 }
00026 
00033 static iBase_EntityHandle embedWithinWorld( bool positive, iGeom_Instance& igm, double world_size, 
00034                                             iBase_EntityHandle body, bool bound_with_world )
00035 {
00036   iBase_EntityHandle final_body;
00037   int igm_result;
00038 
00039   if( !positive && !bound_with_world ){
00040     final_body = body;
00041   }
00042   else{
00043     iBase_EntityHandle world_sphere = makeWorldSphere( igm, world_size );    
00044     
00045     if( positive ){
00046       iGeom_subtractEnts( igm, world_sphere, body, &final_body, &igm_result);
00047       CHECK_IGEOM( igm_result, "making positive body" );
00048     }
00049     else{ // !positive && bound_with_world
00050       iGeom_intersectEnts( igm, world_sphere, body, &final_body, &igm_result);
00051       CHECK_IGEOM( igm_result, "making negative body" );
00052     }
00053   }
00054   return final_body;
00055 }
00056 
00057 
00058 
00059 iBase_EntityHandle applyTransform( const Transform& t, iGeom_Instance& igm, iBase_EntityHandle& e ) {
00060   
00061   int igm_result;
00062   if( t.hasRot()  ){
00063     const Vector3d& axis = t.getAxis();
00064     iGeom_rotateEnt( igm, e, t.getTheta(), axis.v[0], axis.v[1], axis.v[2], &igm_result );
00065     CHECK_IGEOM( igm_result, "applying rotation" );
00066   }
00067   
00068   if( t.hasInversion() ){
00069       //TODO: ask about the right way to do this.  Rotate seems to work, but I don't really know why...
00070       //iGeom_rotateEnt( igm, e, 180, 0, 0, 0, &igm_result );
00071 
00072     //if( !t.hasRot() ){
00073       iGeom_reflectEnt( igm, e, 0, 0, 0, 0, 0, 1, &igm_result );
00074       iGeom_reflectEnt( igm, e, 0, 0, 0, 0, 1, 0, &igm_result );
00075       iGeom_reflectEnt( igm, e, 0, 0, 0, 1, 0, 0, &igm_result );
00076       //}
00077       //else{
00078       //   const Vector3d& axis = t.getAxis();
00079       //   iGeom_reflectEnt( igm, e, axis.v[0], axis.v[1], axis.v[2], &igm_result );
00080       //   }
00081 
00082       CHECK_IGEOM( igm_result, "inverting for transformation" );
00083     }
00084 
00085 
00086   const Vector3d& translation = t.getTranslation();
00087   iGeom_moveEnt( igm, e, translation.v[0], translation.v[1], translation.v[2], &igm_result);
00088   CHECK_IGEOM( igm_result, "applying translation" );
00089   
00090   return e;
00091 }
00092 
00093 iBase_EntityHandle applyReverseTransform( const Transform& tx, iGeom_Instance& igm, iBase_EntityHandle& e ) {
00094   
00095   int igm_result;
00096   Transform rev_t = tx.reverse();
00097 
00098   const Vector3d& translation = rev_t.getTranslation();
00099   iGeom_moveEnt( igm, e, translation.v[0], translation.v[1], translation.v[2], &igm_result);
00100   CHECK_IGEOM( igm_result, "applying reverse translation" );
00101 
00102   if( rev_t.hasInversion() ){
00103     iGeom_rotateEnt( igm, e, 180, 0, 0, 0, &igm_result );
00104     CHECK_IGEOM( igm_result, "inverting for reverse transformation" );
00105   }
00106 
00107   if( rev_t.hasRot() ){
00108     const Vector3d& axis = rev_t.getAxis();
00109     iGeom_rotateEnt( igm, e, rev_t.getTheta(), axis.v[0], axis.v[1], axis.v[2], &igm_result );
00110     CHECK_IGEOM( igm_result, "applying rotation" );
00111   }
00112   
00113   
00114   return e;
00115 }
00116 
00117 
00118 iBase_EntityHandle SurfaceVolume::define( bool positive, iGeom_Instance& igm, double world_size ){
00119   iBase_EntityHandle handle = this->getHandle( positive, igm, world_size );
00120   if( transform ){
00121     handle = applyTransform( *transform, igm, handle );
00122   }
00123   return handle;
00124 }
00125 
00126 
00127 
00128 class PlaneSurface : public SurfaceVolume { 
00129 
00130 protected:
00131   Vector3d normal;
00132   double offset;
00133 
00134 public:
00135   PlaneSurface( const Vector3d& normal_p, double offset_p ) :
00136     SurfaceVolume(), normal(normal_p), offset(offset_p) 
00137   {}
00138   
00139   virtual double getFarthestExtentFromOrigin() const{
00140     // this is a funny situation, since planes are technically infinte...
00141     // in order to have a sane answer, we just return the offset from the origin.
00142     // (multiplied by root 3, which was done in the old converter, why?)
00143     return sqrt(3.0) *  std::fabs(offset);
00144   }
00145 
00146 protected:
00147   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size){
00148 
00149     int igm_result;
00150     iBase_EntityHandle world_sphere = makeWorldSphere(igm, world_size);
00151     iBase_EntityHandle hemisphere;
00152     // note the reversal of sense in this call; mcnp and igeom define it differently.
00153     iGeom_sectionEnt( igm, world_sphere, 
00154                       normal.v[0], normal.v[1], normal.v[2], offset, !positive, &hemisphere, &igm_result);
00155     CHECK_IGEOM( igm_result, "Sectioning world for a plane" );
00156     return hemisphere;
00157 
00158 
00159   }
00160 
00161 };
00162 
00163 
00164 typedef  enum { X=0, Y=1, Z=2 } axis_t;
00165 
00166 
00167 class CylinderSurface : public SurfaceVolume {
00168 
00169 protected:
00170   axis_t axis;
00171   double radius;
00172   Vector3d center;
00173   bool onaxis;
00174 
00175 public:
00176   CylinderSurface( axis_t axis_p, double radius_p ):
00177     SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(true)
00178   {}
00179 
00180   CylinderSurface( axis_t axis_p, double radius_p, double trans1, double trans2 ):
00181     SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(false)
00182   {
00183     switch(axis){
00184     case X: center.v[Y] += trans1; center.v[Z] += trans2; break;
00185     case Y: center.v[X] += trans1; center.v[Z] += trans2; break;
00186     case Z: center.v[X] += trans1; center.v[Y] += trans2; break;
00187     }
00188   }
00189   
00190   virtual double getFarthestExtentFromOrigin( ) const{
00191     return radius + center.length();
00192   }
00193 
00194 protected:
00195   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00196     int igm_result;
00197 
00198     iBase_EntityHandle cylinder;
00199     iGeom_createCylinder( igm, 2.0 * world_size, radius, 0, &cylinder, &igm_result);
00200     CHECK_IGEOM( igm_result, "making cylinder" );
00201 
00202     
00203     if( axis == X ){
00204       iGeom_rotateEnt( igm, cylinder, 90, 0, 1, 0, &igm_result );
00205       CHECK_IGEOM( igm_result, "rotating cylinder (X)" );
00206     }
00207     else if( axis == Y ){
00208       iGeom_rotateEnt( igm, cylinder, 90, 1, 0, 0, &igm_result );
00209       CHECK_IGEOM( igm_result, "rotating cylinder (Y)" );
00210     }
00211 
00212     if( onaxis == false ){
00213       iGeom_moveEnt( igm, cylinder, center.v[0], center.v[1], center.v[2], &igm_result);
00214       CHECK_IGEOM( igm_result, "moving cylinder" );
00215     }
00216 
00217     iBase_EntityHandle final_cylinder = embedWithinWorld( positive, igm, world_size, cylinder, true );
00218 
00219     return final_cylinder;
00220   };
00221 
00222 };
00223 
00224 //#ifdef HAVE_IGEOM_CONE
00225 
00226 class ConeSurface : public SurfaceVolume {
00227 
00228 protected:
00229 
00230   enum nappe {LEFT=-1, BOTH=0, RIGHT=1};
00231 
00232   static enum nappe make_nappe( double param ){
00233     
00234     enum nappe n = static_cast<enum nappe>(param);
00235     if( -1 <= n && n <= 1 ){
00236       return n;
00237     }
00238     else{
00239       std::cerr << "WARNING: Bad cylinder +/-1 argument: " << param << std::endl;
00240       std::cerr << "         will pretend it was really 0" << std::endl;
00241       return BOTH;
00242     }
00243   }
00244 
00245   axis_t axis;
00246   double theta;    
00247   Vector3d center; 
00248   bool onaxis;
00249   enum nappe nappe; 
00250 
00251 public:
00252   ConeSurface( axis_t axis_p, double tsquared_p, double point_p, double nappe_p ):
00253     SurfaceVolume(), axis(axis_p), theta( atan(sqrt(tsquared_p)) ), center(origin), onaxis(true), nappe(make_nappe(nappe_p))
00254   {
00255     center.v[axis] = point_p; 
00256   }
00257 
00258   ConeSurface( axis_t axis_p, double tsquared_p, Vector3d center_p, double nappe_p ):
00259     SurfaceVolume(), axis(axis_p), theta( atan(sqrt(tsquared_p)) ), center(center_p), onaxis(false), nappe(make_nappe(nappe_p))
00260   {}
00261   
00262   virtual double getFarthestExtentFromOrigin( ) const{
00263     return sqrt(3) + ( center.length() );
00264   }
00265 
00266 protected:
00267   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00268 
00269     double height = (center.length() + world_size);
00270 
00271     // based on the textual descriptions in the manual, I think the following expression should be 
00272     // height * tan ( theta / 2 ) -- unless "opening angle" refers to only half the apex angle
00273     // of the cylinder.  But this implementation seems to be more correct in examples I can check against.
00274     double base_radius = height * tan( theta );
00275 
00276     int igm_result;
00277 
00278     iBase_EntityHandle right_nappe = 0;
00279     iBase_EntityHandle left_nappe = 0;
00280     iBase_EntityHandle cone; 
00281     
00282     if( nappe != LEFT){
00283       iGeom_createCone( igm, height, base_radius, 0, 0, &right_nappe, &igm_result);
00284       CHECK_IGEOM( igm_result, "making cone (right nappe)" );
00285       iGeom_rotateEnt( igm, right_nappe, 180, 1, 0, 0, &igm_result);
00286       CHECK_IGEOM( igm_result, "Rotating cone (right nappe)");
00287       iGeom_moveEnt( igm, right_nappe, 0, 0, height/2.0, &igm_result );
00288       CHECK_IGEOM( igm_result, "Moving cone (right nappe)");      
00289       cone = right_nappe;
00290     }
00291     if( nappe != RIGHT ){
00292       iGeom_createCone( igm, height, base_radius, 0, 0, &left_nappe, &igm_result );
00293       CHECK_IGEOM( igm_result, "making cone (left nappe)" );
00294       iGeom_moveEnt( igm, left_nappe, 0, 0, -height/2.0, &igm_result );
00295       CHECK_IGEOM( igm_result, "Moving cone (left nappe)" );
00296       cone = left_nappe;
00297     }
00298 
00299     
00300     if( right_nappe && left_nappe ){
00301       iBase_EntityHandle nappes[2] = {right_nappe, left_nappe};
00302       iGeom_uniteEnts( igm, nappes, 2, &cone, &igm_result );
00303       CHECK_IGEOM( igm_result, "Unioning cone nappes" );
00304     }
00305 
00306     if( axis == X ){
00307       iGeom_rotateEnt( igm, cone, 90, 0, 1, 0, &igm_result );
00308       CHECK_IGEOM( igm_result, "rotating cone (X)" );
00309     }
00310     else if( axis == Y ){
00311       iGeom_rotateEnt( igm, cone, -90, 1, 0, 0, &igm_result );
00312       CHECK_IGEOM( igm_result, "rotating cone (Y)" );
00313     }
00314 
00315     iGeom_moveEnt( igm, cone, center.v[0], center.v[1], center.v[2], &igm_result);
00316     CHECK_IGEOM( igm_result, "moving cone to its apex" );
00317 
00318     iBase_EntityHandle final_cone = embedWithinWorld( positive, igm, world_size, cone, true );
00319 
00320     return final_cone;
00321 
00322     }
00323 
00324 };
00325 
00326 //#endif /* HAVE_IGEOM_CONE */
00327 
00328 class TorusSurface : public SurfaceVolume {
00329 
00330 protected:
00331   Vector3d center;
00332   axis_t axis;
00333   double radius;
00334   double ellipse_axis_rad;
00335   double ellipse_perp_rad;
00336 
00337 public:
00338   TorusSurface( axis_t axis_p, const Vector3d& center_p, double A, double B, double C ) :
00339     center(center_p), axis( axis_p ), radius(A), ellipse_axis_rad( B ), ellipse_perp_rad( C )
00340   {}
00341 
00342   virtual double getFarthestExtentFromOrigin ( ) const {
00343     return center.length() + radius + std::max( ellipse_axis_rad, ellipse_perp_rad );
00344   }
00345 
00346 protected:
00347   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00348 
00349     int igm_result;
00350 
00351     iBase_EntityHandle torus;
00352 
00353     iGeom_createTorus( igm, radius, ellipse_perp_rad, &torus, &igm_result );
00354     CHECK_IGEOM( igm_result, "Creating initial torus");
00355 
00356     if( ellipse_axis_rad != ellipse_perp_rad ){
00357       double scalef = ellipse_axis_rad / ellipse_perp_rad;
00358       iGeom_scaleEnt( igm, torus, 0, 0, 0, 1.0, 1.0, scalef, &igm_result );
00359       CHECK_IGEOM( igm_result, "Scaling torus" );
00360     }
00361     
00362     if( axis == X ){
00363       iGeom_rotateEnt( igm, torus, 90, 0, 1, 0, &igm_result );
00364       CHECK_IGEOM( igm_result, "rotating torus (X)" );
00365     }
00366     else if( axis == Y ){
00367       iGeom_rotateEnt( igm, torus, -90, 1, 0, 0, &igm_result );
00368       CHECK_IGEOM( igm_result, "rotating torus (Y)" );
00369     }
00370 
00371     iGeom_moveEnt( igm, torus, center.v[0], center.v[1], center.v[2], &igm_result);
00372     CHECK_IGEOM( igm_result, "moving torus to its center point" );
00373     
00374     
00375     iBase_EntityHandle final_torus = embedWithinWorld( positive, igm, world_size, torus, false );
00376 
00377     return final_torus;
00378 
00379   }
00380 
00381 };
00382 
00383 class SphereSurface : public SurfaceVolume {
00384 
00385 protected:
00386   Vector3d center;
00387   double radius;
00388 
00389 public:
00390   SphereSurface( const Vector3d& center_p, double radius_p ) :
00391     SurfaceVolume(), center(center_p), radius(radius_p)
00392   {}
00393 
00394   virtual ~SphereSurface(){}
00395 
00396   virtual double getFarthestExtentFromOrigin ( ) const {
00397     return (center.length() + radius);
00398   }
00399 
00400 protected:
00401   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00402 
00403     int igm_result;
00404     iBase_EntityHandle sphere;
00405 
00406     iGeom_createSphere( igm, radius, &sphere, &igm_result);
00407     CHECK_IGEOM( igm_result, "making sphere" );
00408 
00409     iGeom_moveEnt( igm, sphere, center.v[0], center.v[1], center.v[2], &igm_result );
00410     CHECK_IGEOM( igm_result, "moving sphere" );
00411 
00412 
00413     iBase_EntityHandle final_sphere = embedWithinWorld( positive, igm, world_size, sphere, false );
00414     
00415     return final_sphere; 
00416   }
00417 
00418 };
00419 
00420 class EllipsoidSurface : public SurfaceVolume {
00421 
00422 protected:
00423   Vector3d center;
00424   Vector3d axes;
00425 
00426 public:
00427   EllipsoidSurface( const Vector3d& center_p, const Vector3d& axes_p ) :
00428     SurfaceVolume(), center(center_p), axes(axes_p)
00429   {}
00430 
00431   virtual ~EllipsoidSurface(){}
00432 
00433   virtual double getFarthestExtentFromOrigin ( ) const {
00434     return (center.length() + axes.length());
00435   }
00436 
00437 protected:
00438   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00439 
00440     int igm_result;
00441     iBase_EntityHandle sphere;
00442     double radius = 1;
00443 
00444     iGeom_createSphere( igm, radius, &sphere, &igm_result);
00445     CHECK_IGEOM( igm_result, "making sphere" );
00446 
00447     iGeom_scaleEnt( igm, sphere, 0, 0, 0, sqrt(1/axes.v[0]), sqrt(1/axes.v[1]), sqrt(1/axes.v[2]), &igm_result);
00448     CHECK_IGEOM( igm_result, "scaling sphere to ellipsoid" );
00449 
00450     iGeom_moveEnt( igm, sphere, center.v[0], center.v[1], center.v[2], &igm_result );
00451     CHECK_IGEOM( igm_result, "moving sphere" );
00452 
00453 
00454     iBase_EntityHandle final_sphere = embedWithinWorld( positive, igm, world_size, sphere, false );
00455     
00456     return final_sphere; 
00457   }
00458 
00459 };
00460 
00461 
00462 
00463 static Transform axesImage( const Vector3d& v1, const Vector3d& v2, const Vector3d &v3, const Vector3d& translation = Vector3d() )
00464 {
00465   Vector3d x(1, 0, 0), y(0, 1, 0), z(0, 0, 1);
00466   Vector3d a1 = v1.normalize(), a2 = v2.normalize(), a3 = v3.normalize();
00467     
00468   if( OPT_DEBUG ) std::cout << "Axes image: " << a1 << " : " << a2 << " : " << a3 << std::endl;
00469 
00470   double rot_matrix[9] = 
00471     { a1.dot(x), a2.dot(x), a3.dot(x),
00472       a1.dot(y), a2.dot(y), a3.dot(y),
00473       a1.dot(z), a2.dot(z), a3.dot(z) };
00474   
00475 
00476   return Transform( rot_matrix, translation );
00477 }
00478 
00479 static Transform imageZAxisTo( const Vector3d& v, const Vector3d& translation = Vector3d() ){
00480   // approach: find two vectors perpendicular to v, then use axesImage.
00481   Vector3d a = v.normalize();
00482 
00483   Vector3d x(1, 0, 0);
00484   Vector3d b = x.cross( v.normalize() );
00485   
00486   if( b.length() < DBL_EPSILON ){
00487     // v is indistinguishable from the x axis
00488     b = Vector3d(0, -1, 0);
00489     if( OPT_DEBUG ) std::cout << "imageZAxisTo: Changing v " << std::endl;
00490   }
00491 
00492   Vector3d c = b.cross( v.normalize() );
00493   return axesImage( c, b, a, translation );
00494   
00495 
00496 }
00497 
00498 
00499 
00500 class BoxVolume : public SurfaceVolume {
00501 
00502 protected:
00503   Vector3d dimensions;
00504   Transform transform;
00505 
00506 public:
00507   BoxVolume( const Vector3d& corner, const Vector3d& v1, const Vector3d& v2, const Vector3d& v3 ) :
00508     dimensions( v1.length(), v2.length(), v3.length() ), transform( axesImage(v1,v2,v3,corner) )
00509   {}
00510   
00511   virtual double getFarthestExtentFromOrigin ( ) const {
00512     return transform.getTranslation().length() + dimensions.length();
00513   }
00514   
00515 protected:
00516   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00517 
00518 
00519     int igm_result;
00520     iBase_EntityHandle box;
00521 
00522     iGeom_createBrick( igm, dimensions.v[0], dimensions.v[1], dimensions.v[2], &box, &igm_result );
00523     CHECK_IGEOM( igm_result, "making box" );
00524 
00525     Vector3d halfdim = dimensions.scale( 1.0 / 2.0 );
00526     iGeom_moveEnt( igm, box, halfdim.v[0], halfdim.v[1], halfdim.v[2], &igm_result );
00527     CHECK_IGEOM( igm_result, "moving box (halfdim)" );
00528 
00529     box = applyTransform( transform, igm, box );
00530     iBase_EntityHandle final_box = embedWithinWorld( positive, igm, world_size, box, false );
00531 
00532     return final_box;
00533   }
00534 
00535 };
00536 
00537 
00538 class RppVolume : public SurfaceVolume {
00539 
00540 protected:
00541   Vector3d dimensions;
00542   Vector3d center_offset;
00543 
00544 public:
00545   RppVolume( const Vector3d& lower, const Vector3d& upper ) 
00546   {
00547     for( int i = 0; i < 3; ++i ){
00548       dimensions.v[i]    = upper.v[i] - lower.v[i];
00549       center_offset.v[i] = (upper.v[i]+lower.v[i]) / 2.0;
00550     }
00551   }
00552   
00553   virtual double getFarthestExtentFromOrigin ( ) const {
00554     return dimensions.length() + center_offset.length();
00555   }
00556   
00557 protected:
00558   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00559 
00560     int igm_result;
00561     iBase_EntityHandle rpp;
00562 
00563     iGeom_createBrick( igm, dimensions.v[0], dimensions.v[1], dimensions.v[2], &rpp, &igm_result );
00564     CHECK_IGEOM( igm_result, "making rpp" );
00565 
00566     iGeom_moveEnt( igm, rpp, center_offset.v[0], center_offset.v[1], center_offset.v[2], &igm_result );
00567     CHECK_IGEOM( igm_result, "moving rpp" );
00568 
00569     iBase_EntityHandle final_rpp = embedWithinWorld( positive, igm, world_size, rpp, false );
00570     return final_rpp;
00571   }
00572 
00573 };
00574 
00575 
00576 class RecVolume : public SurfaceVolume { 
00577 
00578 protected:
00579   Vector3d base_center;
00580   Transform transform;
00581   double length, radius1, radius2;
00582 
00583 public:
00584   RecVolume( const Vector3d& center_p, const Vector3d& axis, const Vector3d& v1, const Vector3d& v2 ) :
00585     base_center( center_p ), transform( axesImage( v1, v2, axis, center_p ) ), 
00586     length( axis.length() ), radius1( v1.length() ), radius2( v2.length() )
00587   {}
00588 
00589   RecVolume( const Vector3d& center_p, const Vector3d& axis, const Vector3d& v1, double length2 ) :
00590     base_center( center_p ), transform( axesImage( v1, v1.cross(axis), axis, center_p ) ), 
00591     length( axis.length() ), radius1( v1.length() ), radius2( length2 )
00592   {}
00593 
00594   virtual double getFarthestExtentFromOrigin ( ) const {
00595     return base_center.length() + length + std::max( radius1, radius2 );
00596   }
00597 
00598 protected:
00599   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00600     int igm_result;
00601     iBase_EntityHandle rec;
00602     
00603     iGeom_createCylinder( igm, length, radius1, radius2, &rec, &igm_result );
00604     CHECK_IGEOM( igm_result, "creating rec" );
00605     
00606 
00607     double movement_factor = length / 2.0;
00608     iGeom_moveEnt( igm, rec, 0, 0, movement_factor, &igm_result );
00609     CHECK_IGEOM( igm_result, "moving rec" );
00610     
00611     rec = applyTransform( transform, igm, rec );
00612     iBase_EntityHandle final_rec = embedWithinWorld( positive, igm, world_size, rec, false );
00613     return final_rec;
00614   }
00615 };
00616 
00617 
00618 class RccVolume : public SurfaceVolume {
00619 
00620 protected:
00621   Vector3d base_center;
00622   Transform transform;
00623   double length, radius;
00624 
00625 public:
00626   RccVolume( const Vector3d& center_p, const Vector3d& axis, double radius_p ) :
00627     base_center( center_p ), transform( imageZAxisTo( axis, center_p ) ), length( axis.length() ), radius(radius_p) 
00628   {}
00629 
00630   virtual double getFarthestExtentFromOrigin ( ) const {
00631     return base_center.length() + length + radius;
00632   }
00633 
00634 protected:
00635   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00636     int igm_result;
00637     iBase_EntityHandle rcc;
00638 
00639     iGeom_createCylinder( igm, length, radius, 0, &rcc, &igm_result );
00640     CHECK_IGEOM( igm_result, "creating rcc" );
00641     
00642     double movement_factor = length / 2.0;
00643     iGeom_moveEnt( igm, rcc, 0, 0, movement_factor, &igm_result );
00644     CHECK_IGEOM( igm_result, "moving rcc" );
00645 
00646     rcc = applyTransform( transform, igm, rcc );
00647     iBase_EntityHandle final_rcc = embedWithinWorld( positive, igm, world_size, rcc, false );
00648     return final_rcc;
00649   }
00650 
00651 };
00652 
00653 
00654 #ifdef HAVE_IGEOM_CONE
00655 class TrcVolume : public SurfaceVolume { 
00656 
00657 protected:
00658   Vector3d base_center;
00659   Transform transform;
00660   double length, radius1, radius2;
00661 
00662 public:
00663   TrcVolume( const Vector3d& center_p, const Vector3d& axis, double radius1_p, double radius2_p ) :
00664     base_center( center_p ), transform( imageZAxisTo( axis, center_p ) ), length( axis.length() ), 
00665     radius1(radius1_p), radius2(radius2_p) 
00666   {}
00667 
00668   virtual double getFarthestExtentFromOrigin ( ) const {
00669     return base_center.length() + length + std::max(radius1,radius2);
00670   }
00671 
00672 protected:
00673   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00674     int igm_result;
00675     iBase_EntityHandle trc;
00676 
00677     iGeom_createCone( igm, length, radius1, 0, radius2, &trc, &igm_result );
00678     CHECK_IGEOM( igm_result, "creating trc" );
00679     
00680     iGeom_moveEnt( igm, trc, 0, 0, length / 2.0, &igm_result );
00681     CHECK_IGEOM( igm_result, "moving trc" );
00682 
00683     trc = applyTransform( transform, igm, trc );
00684     iBase_EntityHandle final_trc = embedWithinWorld( positive, igm, world_size, trc, false );
00685     return final_trc;
00686   }
00687 
00688 };
00689 #endif
00690 
00691 
00692 class HexVolume : public SurfaceVolume { 
00693 
00694 protected:
00695   Vector3d base_center;
00696   Vector3d heightV, RV, SV, TV; 
00697   //Transform transform;
00698 
00699 public:
00700   HexVolume( const Vector3d& center_p, const Vector3d& h_p, const Vector3d& r_p, 
00701              const Vector3d& s_p, const Vector3d& t_p ) :
00702     base_center(center_p), heightV(h_p), RV(r_p), SV(s_p), TV(t_p)
00703   {}
00704 
00705   HexVolume( const Vector3d& center_p, const Vector3d& h_p, const Vector3d& r_p ) :
00706     base_center(center_p), heightV(h_p), RV(r_p), SV( r_p.rotate_about(h_p,60.0) ), TV( r_p.rotate_about(h_p,120.0) )
00707   {
00708     if( OPT_DEBUG ){ 
00709       std::cout << "Inferred vectors for 9-args HEX/RHP:" << RV << SV << TV << std::endl;
00710     }
00711   }
00712 
00713   virtual double getFarthestExtentFromOrigin ( ) const { 
00714     double hex_max = std::max( RV.length(), std::max( SV.length(), TV.length() ) );
00715     return base_center.length() + heightV.length() + hex_max; 
00716   }
00717 
00718 protected:
00719   virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size ){
00720     int igm_result;
00721     iBase_EntityHandle hex;
00722 
00723     hex = makeWorldSphere( igm, world_size );
00724 
00725     Vector3d b = - heightV.normalize();
00726     iGeom_sectionEnt( igm, hex, b.v[0], b.v[1], b.v[2], 0, true, &hex, &igm_result );
00727     CHECK_IGEOM( igm_result, "Sectioning world for a hex (1)" );
00728 
00729     b = -b;
00730     iGeom_sectionEnt( igm, hex, b.v[0], b.v[1], b.v[2], heightV.length(), true, &hex, &igm_result );
00731     CHECK_IGEOM( igm_result, "Sectioning world for a hex (2)" );
00732 
00733 
00734     
00735     const Vector3d* vec[3] = {&RV, &SV, &TV};
00736     for( int i = 0; i < 3; ++i ){
00737       Vector3d v = *(vec[i]);
00738       double length = v.length(); 
00739       v = v.normalize();
00740       
00741       iGeom_sectionEnt( igm, hex, v.v[0], v.v[1], v.v[2], length, true, &hex, &igm_result );
00742       CHECK_IGEOM( igm_result, "Sectioning world for a hex (3)" );
00743       
00744       v = -v;
00745       iGeom_sectionEnt( igm, hex, v.v[0], v.v[1], v.v[2], length, true, &hex, &igm_result );
00746       CHECK_IGEOM( igm_result, "Sectioning world for a hex (4)" );
00747       
00748 
00749     }
00750 
00751     iGeom_moveEnt( igm, hex, base_center.v[0], base_center.v[1], base_center.v[2], &igm_result );
00752     CHECK_IGEOM( igm_result, "Moving hex" );
00753 
00754     iBase_EntityHandle final_hex = embedWithinWorld( positive, igm, world_size, hex, false );
00755     return final_hex;
00756 
00757   }
00758 
00759 };
00760 
00761 
00762 class VolumeCache{
00763 
00764 protected:
00765   std::map<const SurfaceCard*,SurfaceVolume*> mapping;
00766 
00767 public:
00768   VolumeCache(){}
00769   
00770   bool contains( const SurfaceCard* c ) const {
00771     return ( mapping.find(c) != mapping.end() );
00772   }
00773 
00774   SurfaceVolume* get( const SurfaceCard* c ) {
00775     assert( contains( c ) );
00776     return (*(mapping.find(c))).second;
00777   }
00778 
00779   void insert( const SurfaceCard* c, SurfaceVolume* s ){
00780     mapping[c] = s;
00781   }
00782 
00783 };
00784 
00785 bool sqIsEllipsoid(const std::vector< double >& args) {
00786 
00787   bool isEllipsoid = true;
00788   
00789   if ( args.at(3)*args.at(3) + args.at(4)*args.at(4) + args.at(5)*args.at(5) != 0)
00790     isEllipsoid = false;
00791 
00792   if ( args.at(0)*args.at(1)*args.at(2)*args.at(6) >= 0 )
00793     isEllipsoid = false;
00794 
00795   return isEllipsoid;
00796   
00797 }
00798 
00799 // the VolumeCache to use if none is provided to makeSurface() calls.
00800 static VolumeCache default_volume_cache;
00801 
00802 SurfaceVolume& makeSurface( const SurfaceCard* card, VolumeCache* v){
00803   
00804   VolumeCache& cache = default_volume_cache;
00805   if( v != NULL ){
00806     cache = *v;
00807   }
00808 
00809   
00810   if( cache.contains( card ) ){
00811     return *cache.get( card );
00812   }
00813   else{ 
00814     // SurfaceCard variables:  mnemonic, args, coord_xform  
00815     SurfaceVolume* surface;
00816     const std::string& mnemonic = card->getMnemonic();
00817     const std::vector< double >& args = card->getArgs(); 
00818     
00819     if( mnemonic == "so"){
00820       surface = new SphereSurface( origin, args.at(0) );
00821     }
00822     else if( mnemonic == "sx"){
00823       surface = new SphereSurface( Vector3d( args.at(0), 0, 0 ), args.at(1) );
00824     }
00825     else if( mnemonic == "sy"){
00826       surface = new SphereSurface( Vector3d( 0, args.at(0), 0 ), args.at(1) );
00827     }
00828     else if( mnemonic == "sz"){
00829       surface = new SphereSurface( Vector3d( 0, 0, args.at(0) ), args.at(1) );
00830     }
00831     else if( mnemonic == "s" || mnemonic == "sph" ){
00832       surface = new SphereSurface( Vector3d( args ), args.at(3) );
00833     }
00834     else if( mnemonic == "sq" && sqIsEllipsoid(args) ){
00835       surface = new EllipsoidSurface( Vector3d( args.at(7), args.at(8), args.at(9) ), Vector3d( args )* (-1/args.at(6)) );
00836     }
00837     else if( mnemonic == "p"){
00838       if( args.size() == 4 ){
00839         // plane given as ABCD coefficients (Ax+By+Cz-D=0)
00840         Vector3d v(args);
00841         surface = new PlaneSurface( v, args.at(3)/v.length() );
00842       }
00843       else if( args.size() == 9 ){
00844         // plane is given at three points in space
00845         Vector3d v1(args), v2(args,3), v3(args,6);
00846 
00847         // ordering of the points is arbitrary, but make sure v1 is furthest from origin.
00848         if( v1.length() < v2.length() ){
00849           std::swap( v1, v2 );
00850         }
00851         if( v1.length() < v3.length() ){
00852           std::swap( v1, v3 );
00853         }
00854 
00855         // the normal (up to reversal) of the plane 
00856         // The normal may need to be reversed to ensure that the origin
00857         // has negative sense, as required by MCNP
00858         Vector3d normal = v1.add(-v2).cross( v1.add(-v3) ).normalize();
00859 
00860         // If a ray started from the origin and followed the normal vector,
00861         // p is the point where it intersects the plane
00862         Vector3d p = normal.scale( v1.dot(normal) );
00863 
00864         // cos of the angle between v1 and normal 
00865         double angle = normal.dot( v1.normalize() );
00866         
00867         // invert the normal if the angle is > 90 degrees, which indicates
00868         // that reversal is required.
00869         if( angle < 0 ){
00870           normal   = -normal;
00871         }
00872 
00873         //std::cout << normal << " " << p.length() << " : " << angle << std::endl;
00874         surface = new PlaneSurface( normal, p.length() );
00875 
00876 
00877       }
00878       else{ 
00879         throw std::runtime_error( "P surface with unsupported number of params" );
00880       }
00881 
00882     }
00883     else if( mnemonic == "px"){
00884       surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
00885     }
00886     else if( mnemonic == "py"){
00887       surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
00888     }
00889     else if( mnemonic == "pz"){
00890       surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
00891     }
00892     else if( mnemonic == "cx" ){
00893       surface = new CylinderSurface( X, args.at(0) );
00894     }
00895     else if( mnemonic == "cy" ){
00896       surface = new CylinderSurface( Y, args.at(0) );
00897     }
00898     else if( mnemonic == "cz" ){
00899       surface = new CylinderSurface( Z, args.at(0) );
00900     }
00901     else if( mnemonic == "c/x"){
00902       surface = new CylinderSurface( X, args.at(2), args.at(0), args.at(1) );
00903     }
00904     else if( mnemonic == "c/y"){
00905       surface = new CylinderSurface( Y, args.at(2), args.at(0), args.at(1) );
00906     }
00907     else if( mnemonic == "c/z"){
00908       surface = new CylinderSurface( Z, args.at(2), args.at(0), args.at(1) );
00909     }
00910 #ifdef HAVE_IGEOM_CONE
00911     else if( mnemonic == "kx"){
00912       double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
00913       surface = new ConeSurface( X, args.at(1), args.at(0), arg3 );
00914     } 
00915     else if( mnemonic == "ky"){
00916       double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
00917       surface = new ConeSurface( Y, args.at(1), args.at(0), arg3 );
00918     }
00919     else if( mnemonic == "kz"){
00920       double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
00921       surface = new ConeSurface( Z, args.at(1), args.at(0), arg3 );
00922     }
00923     else if( mnemonic == "k/x" ){
00924       double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
00925       surface = new ConeSurface( X, args.at(3), Vector3d(args),  arg5 );
00926     }
00927     else if( mnemonic == "k/y" ){
00928       double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
00929       surface = new ConeSurface( Y, args.at(3), Vector3d(args),  arg5 );
00930     }
00931     else if( mnemonic == "k/z" ){
00932       double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
00933       surface = new ConeSurface( Z, args.at(3), Vector3d(args),  arg5 );
00934     }
00935     else if( mnemonic == "trc" ){
00936       surface = new TrcVolume( Vector3d(args), Vector3d(args,3), args.at(6), args.at(7) );
00937     }
00938 #endif /*HAVE_IGEOM_CONE */
00939     else if( mnemonic == "tx" ){
00940       surface = new TorusSurface( X, Vector3d(args), args.at(3), args.at(4), args.at(5) );
00941     } 
00942     else if( mnemonic == "ty" ){
00943       surface = new TorusSurface( Y, Vector3d(args), args.at(3), args.at(4), args.at(5) );
00944     } 
00945     else if( mnemonic == "tz" ){
00946       surface = new TorusSurface( Z, Vector3d(args), args.at(3), args.at(4), args.at(5) );
00947     }
00948     else if( mnemonic == "box" ){
00949       surface = new BoxVolume( Vector3d(args), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9) );
00950     }
00951     else if( mnemonic == "rpp" ){
00952       surface = new RppVolume( Vector3d(args.at(0), args.at(2), args.at(4)), Vector3d(args.at(1), args.at(3), args.at(5)) );
00953     }
00954     else if( mnemonic == "rcc" ){
00955       surface = new RccVolume( Vector3d(args), Vector3d(args,3), args.at(6) );
00956     }
00957     else if( mnemonic == "rec" ){
00958       if( args.size() == 10 ){
00959         surface = new RecVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6),  args.at(9) );
00960       }
00961       else{
00962         surface = new RecVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9) );
00963       }
00964     }
00965     else if( mnemonic == "hex" || mnemonic == "rhp" ){
00966       if( args.size() == 9 ){
00967         surface = new HexVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6) );
00968       }
00969       else{
00970         surface = new HexVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9), Vector3d(args,12) );
00971       }
00972     }
00973     else if( mnemonic == "x" ) {
00974       switch (args.size()) {
00975       case 2: // plane
00976         surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
00977         break;
00978       case 4: // either plane, cylinder or cone
00979         if ( args.at(0) == args.at(2) ) // plane
00980           surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
00981         else if (args.at(1) == args.at(3)) // cylinder
00982           surface = new CylinderSurface( X, args.at(1) );
00983         else // cone
00984           {
00985             double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
00986             double apex_p = args.at(0) - args.at(1)/m;
00987             surface = new ConeSurface( X, m*m, apex_p, (m > 0 ? 1 : -1 ) );
00988           }
00989         break;
00990       default:
00991         throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
00992         break;
00993       }
00994     }
00995     else if( mnemonic == "y" ) {
00996       switch (args.size()) {
00997       case 2: // plane
00998         surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
00999         break;
01000       case 4: // either plane, cylinder or cone
01001         if ( args.at(0) == args.at(2) ) // plane
01002           surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
01003         else if (args.at(1) == args.at(3)) // cylinder
01004           surface = new CylinderSurface( Y, args.at(1) );
01005         else // cone
01006           {
01007             double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
01008             double apex_p = args.at(0) - args.at(1)/m;
01009             surface = new ConeSurface( Y, m*m, apex_p, (m > 0 ? 1 : -1 ) );
01010           }
01011         break;
01012       default:
01013         throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
01014         break;
01015       }
01016     }
01017     else if( mnemonic == "z" ) {
01018       switch (args.size()) {
01019       case 2: // plane
01020         surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
01021         break;
01022       case 4: // either plane, cylinder or cone
01023         if ( args.at(0) == args.at(2) ) // plane
01024           surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
01025         else if (args.at(1) == args.at(3)) // cylinder
01026           surface = new CylinderSurface( Z, args.at(1) );
01027         else // cone
01028           {
01029             double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
01030             double apex_p = args.at(0) - args.at(1)/m;
01031             surface = new ConeSurface( Z, m*m, apex_p, (m > 0 ? 1 : -1 ) );
01032           }
01033         break;
01034       default:
01035         throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
01036         break;
01037       }
01038     }
01039     else{
01040       throw std::runtime_error( mnemonic + " is not a supported surface" );
01041     }
01042     
01043     if( card->getTransform().hasData() ){
01044       const Transform& transform = card->getTransform().getData();
01045       surface->setTransform( &transform );
01046     }
01047     
01048     cache.insert( card, surface );
01049     return *surface;
01050     
01051   }
01052   
01053   
01054 }
01055 
01056 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines