cgma
volumes.hpp File Reference
#include <cstdlib>
#include "iGeom.h"

Go to the source code of this file.

Classes

class  SurfaceVolume

Defines

#define CHECK_BUF_SIZE   512
#define CHECK_IGEOM(err, msg)

Functions

SurfaceVolumemakeSurface (const SurfaceCard *card, VolumeCache *v=NULL)
iBase_EntityHandle makeWorldSphere (iGeom_Instance &igm, double world_size)
iBase_EntityHandle applyTransform (const Transform &t, iGeom_Instance &igm, iBase_EntityHandle &e)
iBase_EntityHandle applyReverseTransform (const Transform &tx, iGeom_Instance &igm, iBase_EntityHandle &e)

Variables

static char m_buf [CHECK_BUF_SIZE]

Define Documentation

#define CHECK_BUF_SIZE   512

Definition at line 59 of file volumes.hpp.

#define CHECK_IGEOM (   err,
  msg 
)
Value:
do{/*std::cout << msg << std::endl;*/ if((err) != iBase_SUCCESS){     \
    std::cerr << "iGeom error (" << err << "): " << msg << std::endl;   \
    iGeom_getDescription( igm, m_buf, CHECK_BUF_SIZE);                  \
    std::cerr << " * " << m_buf << std::endl;                           \
     } } while(0)

Definition at line 62 of file volumes.hpp.


Function Documentation

Definition at line 93 of file volumes.cpp.

                                                                                                            {
  
  int igm_result;
  Transform rev_t = tx.reverse();

  const Vector3d& translation = rev_t.getTranslation();
  iGeom_moveEnt( igm, e, translation.v[0], translation.v[1], translation.v[2], &igm_result);
  CHECK_IGEOM( igm_result, "applying reverse translation" );

  if( rev_t.hasInversion() ){
    iGeom_rotateEnt( igm, e, 180, 0, 0, 0, &igm_result );
    CHECK_IGEOM( igm_result, "inverting for reverse transformation" );
  }

  if( rev_t.hasRot() ){
    const Vector3d& axis = rev_t.getAxis();
    iGeom_rotateEnt( igm, e, rev_t.getTheta(), axis.v[0], axis.v[1], axis.v[2], &igm_result );
    CHECK_IGEOM( igm_result, "applying rotation" );
  }
  
  
  return e;
}

Definition at line 59 of file volumes.cpp.

                                                                                                    {
  
  int igm_result;
  if( t.hasRot()  ){
    const Vector3d& axis = t.getAxis();
    iGeom_rotateEnt( igm, e, t.getTheta(), axis.v[0], axis.v[1], axis.v[2], &igm_result );
    CHECK_IGEOM( igm_result, "applying rotation" );
  }
  
  if( t.hasInversion() ){
      //TODO: ask about the right way to do this.  Rotate seems to work, but I don't really know why...
      //iGeom_rotateEnt( igm, e, 180, 0, 0, 0, &igm_result );

    //if( !t.hasRot() ){
      iGeom_reflectEnt( igm, e, 0, 0, 0, 0, 0, 1, &igm_result );
      iGeom_reflectEnt( igm, e, 0, 0, 0, 0, 1, 0, &igm_result );
      iGeom_reflectEnt( igm, e, 0, 0, 0, 1, 0, 0, &igm_result );
      //}
      //else{
      //   const Vector3d& axis = t.getAxis();
      //   iGeom_reflectEnt( igm, e, axis.v[0], axis.v[1], axis.v[2], &igm_result );
      //   }

      CHECK_IGEOM( igm_result, "inverting for transformation" );
    }


  const Vector3d& translation = t.getTranslation();
  iGeom_moveEnt( igm, e, translation.v[0], translation.v[1], translation.v[2], &igm_result);
  CHECK_IGEOM( igm_result, "applying translation" );
  
  return e;
}
SurfaceVolume& makeSurface ( const SurfaceCard card,
VolumeCache v = NULL 
)

Function to create an SurfaceVolume object from a SurfaceCard. Created volumes are kept in a cache. If the v parameter is null, a default cache (static within volumes.cpp) will be used. If multiple MCNPInput instances are used in the course of a single program, then there will need to be a way to get multiple VolumeCache objects; that's not an issue right now, so I haven't coded it.

Definition at line 802 of file volumes.cpp.

                                                                    {
  
  VolumeCache& cache = default_volume_cache;
  if( v != NULL ){
    cache = *v;
  }

  
  if( cache.contains( card ) ){
    return *cache.get( card );
  }
  else{ 
    // SurfaceCard variables:  mnemonic, args, coord_xform  
    SurfaceVolume* surface;
    const std::string& mnemonic = card->getMnemonic();
    const std::vector< double >& args = card->getArgs(); 
    
    if( mnemonic == "so"){
      surface = new SphereSurface( origin, args.at(0) );
    }
    else if( mnemonic == "sx"){
      surface = new SphereSurface( Vector3d( args.at(0), 0, 0 ), args.at(1) );
    }
    else if( mnemonic == "sy"){
      surface = new SphereSurface( Vector3d( 0, args.at(0), 0 ), args.at(1) );
    }
    else if( mnemonic == "sz"){
      surface = new SphereSurface( Vector3d( 0, 0, args.at(0) ), args.at(1) );
    }
    else if( mnemonic == "s" || mnemonic == "sph" ){
      surface = new SphereSurface( Vector3d( args ), args.at(3) );
    }
    else if( mnemonic == "sq" && sqIsEllipsoid(args) ){
      surface = new EllipsoidSurface( Vector3d( args.at(7), args.at(8), args.at(9) ), Vector3d( args )* (-1/args.at(6)) );
    }
    else if( mnemonic == "p"){
      if( args.size() == 4 ){
        // plane given as ABCD coefficients (Ax+By+Cz-D=0)
        Vector3d v(args);
        surface = new PlaneSurface( v, args.at(3)/v.length() );
      }
      else if( args.size() == 9 ){
        // plane is given at three points in space
        Vector3d v1(args), v2(args,3), v3(args,6);

        // ordering of the points is arbitrary, but make sure v1 is furthest from origin.
        if( v1.length() < v2.length() ){
          std::swap( v1, v2 );
        }
        if( v1.length() < v3.length() ){
          std::swap( v1, v3 );
        }

        // the normal (up to reversal) of the plane 
        // The normal may need to be reversed to ensure that the origin
        // has negative sense, as required by MCNP
        Vector3d normal = v1.add(-v2).cross( v1.add(-v3) ).normalize();

        // If a ray started from the origin and followed the normal vector,
        // p is the point where it intersects the plane
        Vector3d p = normal.scale( v1.dot(normal) );

        // cos of the angle between v1 and normal 
        double angle = normal.dot( v1.normalize() );
        
        // invert the normal if the angle is > 90 degrees, which indicates
        // that reversal is required.
        if( angle < 0 ){
          normal   = -normal;
        }

        //std::cout << normal << " " << p.length() << " : " << angle << std::endl;
        surface = new PlaneSurface( normal, p.length() );


      }
      else{ 
        throw std::runtime_error( "P surface with unsupported number of params" );
      }

    }
    else if( mnemonic == "px"){
      surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
    }
    else if( mnemonic == "py"){
      surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
    }
    else if( mnemonic == "pz"){
      surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
    }
    else if( mnemonic == "cx" ){
      surface = new CylinderSurface( X, args.at(0) );
    }
    else if( mnemonic == "cy" ){
      surface = new CylinderSurface( Y, args.at(0) );
    }
    else if( mnemonic == "cz" ){
      surface = new CylinderSurface( Z, args.at(0) );
    }
    else if( mnemonic == "c/x"){
      surface = new CylinderSurface( X, args.at(2), args.at(0), args.at(1) );
    }
    else if( mnemonic == "c/y"){
      surface = new CylinderSurface( Y, args.at(2), args.at(0), args.at(1) );
    }
    else if( mnemonic == "c/z"){
      surface = new CylinderSurface( Z, args.at(2), args.at(0), args.at(1) );
    }
#ifdef HAVE_IGEOM_CONE
    else if( mnemonic == "kx"){
      double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
      surface = new ConeSurface( X, args.at(1), args.at(0), arg3 );
    } 
    else if( mnemonic == "ky"){
      double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
      surface = new ConeSurface( Y, args.at(1), args.at(0), arg3 );
    }
    else if( mnemonic == "kz"){
      double arg3 = ( args.size() == 3 ? args.at(2) : 0.0 );
      surface = new ConeSurface( Z, args.at(1), args.at(0), arg3 );
    }
    else if( mnemonic == "k/x" ){
      double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
      surface = new ConeSurface( X, args.at(3), Vector3d(args),  arg5 );
    }
    else if( mnemonic == "k/y" ){
      double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
      surface = new ConeSurface( Y, args.at(3), Vector3d(args),  arg5 );
    }
    else if( mnemonic == "k/z" ){
      double arg5 = ( args.size() == 5 ? args.at(4) : 0.0 );
      surface = new ConeSurface( Z, args.at(3), Vector3d(args),  arg5 );
    }
    else if( mnemonic == "trc" ){
      surface = new TrcVolume( Vector3d(args), Vector3d(args,3), args.at(6), args.at(7) );
    }
#endif /*HAVE_IGEOM_CONE */
    else if( mnemonic == "tx" ){
      surface = new TorusSurface( X, Vector3d(args), args.at(3), args.at(4), args.at(5) );
    } 
    else if( mnemonic == "ty" ){
      surface = new TorusSurface( Y, Vector3d(args), args.at(3), args.at(4), args.at(5) );
    } 
    else if( mnemonic == "tz" ){
      surface = new TorusSurface( Z, Vector3d(args), args.at(3), args.at(4), args.at(5) );
    }
    else if( mnemonic == "box" ){
      surface = new BoxVolume( Vector3d(args), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9) );
    }
    else if( mnemonic == "rpp" ){
      surface = new RppVolume( Vector3d(args.at(0), args.at(2), args.at(4)), Vector3d(args.at(1), args.at(3), args.at(5)) );
    }
    else if( mnemonic == "rcc" ){
      surface = new RccVolume( Vector3d(args), Vector3d(args,3), args.at(6) );
    }
    else if( mnemonic == "rec" ){
      if( args.size() == 10 ){
        surface = new RecVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6),  args.at(9) );
      }
      else{
        surface = new RecVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9) );
      }
    }
    else if( mnemonic == "hex" || mnemonic == "rhp" ){
      if( args.size() == 9 ){
        surface = new HexVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6) );
      }
      else{
        surface = new HexVolume( Vector3d( args ), Vector3d(args,3), Vector3d(args,6), Vector3d(args,9), Vector3d(args,12) );
      }
    }
    else if( mnemonic == "x" ) {
      switch (args.size()) {
      case 2: // plane
        surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
        break;
      case 4: // either plane, cylinder or cone
        if ( args.at(0) == args.at(2) ) // plane
          surface = new PlaneSurface( Vector3d( 1, 0, 0), args.at(0) );
        else if (args.at(1) == args.at(3)) // cylinder
          surface = new CylinderSurface( X, args.at(1) );
        else // cone
          {
            double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
            double apex_p = args.at(0) - args.at(1)/m;
            surface = new ConeSurface( X, m*m, apex_p, (m > 0 ? 1 : -1 ) );
          }
        break;
      default:
        throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
        break;
      }
    }
    else if( mnemonic == "y" ) {
      switch (args.size()) {
      case 2: // plane
        surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
        break;
      case 4: // either plane, cylinder or cone
        if ( args.at(0) == args.at(2) ) // plane
          surface = new PlaneSurface( Vector3d( 0, 1, 0), args.at(0) );
        else if (args.at(1) == args.at(3)) // cylinder
          surface = new CylinderSurface( Y, args.at(1) );
        else // cone
          {
            double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
            double apex_p = args.at(0) - args.at(1)/m;
            surface = new ConeSurface( Y, m*m, apex_p, (m > 0 ? 1 : -1 ) );
          }
        break;
      default:
        throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
        break;
      }
    }
    else if( mnemonic == "z" ) {
      switch (args.size()) {
      case 2: // plane
        surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
        break;
      case 4: // either plane, cylinder or cone
        if ( args.at(0) == args.at(2) ) // plane
          surface = new PlaneSurface( Vector3d( 0, 0, 1), args.at(0) );
        else if (args.at(1) == args.at(3)) // cylinder
          surface = new CylinderSurface( Z, args.at(1) );
        else // cone
          {
            double m = (args.at(3) - args.at(1))/(args.at(2)-args.at(0));
            double apex_p = args.at(0) - args.at(1)/m;
            surface = new ConeSurface( Z, m*m, apex_p, (m > 0 ? 1 : -1 ) );
          }
        break;
      default:
        throw std::runtime_error( mnemonic + " is only a supported surface with 2 or 4 arguments" );
        break;
      }
    }
    else{
      throw std::runtime_error( mnemonic + " is not a supported surface" );
    }
    
    if( card->getTransform().hasData() ){
      const Transform& transform = card->getTransform().getData();
      surface->setTransform( &transform );
    }
    
    cache.insert( card, surface );
    return *surface;
    
  }
  
  
}
iBase_EntityHandle makeWorldSphere ( iGeom_Instance igm,
double  world_size 
)

Definition at line 16 of file volumes.cpp.

                                                                            {
  iBase_EntityHandle world_sphere;
  int igm_result;
  // Note: I tried using createBrick instead of createSphere to bound the universe with a box
  // instead of a sphere.  This worked but led to a substantial increase in run times and
  // memory usage, so should be avoided.
  iGeom_createSphere( igm, world_size, &world_sphere, &igm_result);
  CHECK_IGEOM( igm_result, "making world sphere" );
  return world_sphere;
}

Variable Documentation

char m_buf[CHECK_BUF_SIZE] [static]

Definition at line 60 of file volumes.hpp.

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines