<commit_before>#include <iostream>
#include <map>
#include <cfloat>

#include <cassert>

#include "MCNPInput.hpp"
#include "volumes.hpp"
#include "geometry.hpp"
#include "options.hpp"


static Vector3d origin(0,0,0);


iBase_EntityHandle makeWorldSphere( iGeom_Instance& igm, double world_size ){
  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;
}

/**
 * A convenience function for SurfaceVolumes to call at the end of getHandle functions.
 * Return the negative or positive sense of the given body, as appropriate.  If
 * bound_within_world is true, a negative-sense body will be intersected with the world sphere
 * (a step necessary for cylinders and other infinite bodies) 
 */
static iBase_EntityHandle embedWithinWorld( bool positive, iGeom_Instance& igm, double world_size, 
                                            iBase_EntityHandle body, bool bound_with_world )
{
  iBase_EntityHandle final_body;
  int igm_result;

  if( !positive && !bound_with_world ){
    final_body = body;
  }
  else{
    iBase_EntityHandle world_sphere = makeWorldSphere( igm, world_size );    
    
    if( positive ){
      iGeom_subtractEnts( igm, world_sphere, body, &final_body, &igm_result);
      CHECK_IGEOM( igm_result, "making positive body" );
    }
    else{ // !positive && bound_with_world
      iGeom_intersectEnts( igm, world_sphere, body, &final_body, &igm_result);
      CHECK_IGEOM( igm_result, "making negative body" );
    }
  }
  return final_body;
}



iBase_EntityHandle applyTransform( const Transform& t, iGeom_Instance& igm, iBase_EntityHandle& e ) {
  
  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;
}

iBase_EntityHandle applyReverseTransform( const Transform& tx, iGeom_Instance& igm, iBase_EntityHandle& e ) {
  
  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;
}


iBase_EntityHandle SurfaceVolume::define( bool positive, iGeom_Instance& igm, double world_size ){
  iBase_EntityHandle handle = this->getHandle( positive, igm, world_size );
  if( transform ){
    handle = applyTransform( *transform, igm, handle );
  }
  return handle;
}



class PlaneSurface : public SurfaceVolume { 

protected:
  Vector3d normal;
  double offset;

public:
  PlaneSurface( const Vector3d& normal_p, double offset_p ) :
    SurfaceVolume(), normal(normal_p), offset(offset_p) 
  {}
  
  virtual double getFarthestExtentFromOrigin() const{
    // this is a funny situation, since planes are technically infinte...
    // in order to have a sane answer, we just return the offset from the origin.
    // (multiplied by root 3, which was done in the old converter, why?)
    return sqrt(3.0) *  std::fabs(offset);
  }

protected:
  virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size){

    int igm_result;
    iBase_EntityHandle world_sphere = makeWorldSphere(igm, world_size);
    iBase_EntityHandle hemisphere;
    // note the reversal of sense in this call; mcnp and igeom define it differently.
    iGeom_sectionEnt( igm, world_sphere, 
                      normal.v[0], normal.v[1], normal.v[2], offset, !positive, &hemisphere, &igm_result);
    CHECK_IGEOM( igm_result, "Sectioning world for a plane" );
    return hemisphere;


  }

};


typedef  enum { X=0, Y=1, Z=2 } axis_t;


class CylinderSurface : public SurfaceVolume {

protected:
  axis_t axis;
  double radius;
  Vector3d center;
  bool onaxis;

public:
  CylinderSurface( axis_t axis_p, double radius_p ):
    SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(true)
  {}

  CylinderSurface( axis_t axis_p, double radius_p, double trans1, double trans2 ):
    SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(false)
  {
    switch(axis){
    case X: center.v[Y] += trans1; center.v[Z] += trans2; break;
    case Y: center.v[X] += trans1; center.v[Z] += trans2; break;
    case Z: center.v[X] += trans1; center.v[Y] += trans2; break;
    }
  }
  
  virtual double getFarthestExtentFromOrigin( ) const{
    return radius + center.length();<commit_msg>Adding the gq class.<commit_after>#include <iostream>
#include <map>
#include <cfloat>

#include <cassert>

#include "MCNPInput.hpp"
#include "volumes.hpp"
#include "geometry.hpp"
#include "options.hpp"


static Vector3d origin(0,0,0);


iBase_EntityHandle makeWorldSphere( iGeom_Instance& igm, double world_size ){
  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;
}

/**
 * A convenience function for SurfaceVolumes to call at the end of getHandle functions.
 * Return the negative or positive sense of the given body, as appropriate.  If
 * bound_within_world is true, a negative-sense body will be intersected with the world sphere
 * (a step necessary for cylinders and other infinite bodies) 
 */
static iBase_EntityHandle embedWithinWorld( bool positive, iGeom_Instance& igm, double world_size, 
                                            iBase_EntityHandle body, bool bound_with_world )
{
  iBase_EntityHandle final_body;
  int igm_result;

  if( !positive && !bound_with_world ){
    final_body = body;
  }
  else{
    iBase_EntityHandle world_sphere = makeWorldSphere( igm, world_size );    
    
    if( positive ){
      iGeom_subtractEnts( igm, world_sphere, body, &final_body, &igm_result);
      CHECK_IGEOM( igm_result, "making positive body" );
    }
    else{ // !positive && bound_with_world
      iGeom_intersectEnts( igm, world_sphere, body, &final_body, &igm_result);
      CHECK_IGEOM( igm_result, "making negative body" );
    }
  }
  return final_body;
}



iBase_EntityHandle applyTransform( const Transform& t, iGeom_Instance& igm, iBase_EntityHandle& e ) {
  
  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;
}

iBase_EntityHandle applyReverseTransform( const Transform& tx, iGeom_Instance& igm, iBase_EntityHandle& e ) {
  
  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;
}


iBase_EntityHandle SurfaceVolume::define( bool positive, iGeom_Instance& igm, double world_size ){
  iBase_EntityHandle handle = this->getHandle( positive, igm, world_size );
  if( transform ){
    handle = applyTransform( *transform, igm, handle );
  }
  return handle;
}



class PlaneSurface : public SurfaceVolume { 

protected:
  Vector3d normal;
  double offset;

public:
  PlaneSurface( const Vector3d& normal_p, double offset_p ) :
    SurfaceVolume(), normal(normal_p), offset(offset_p) 
  {}
  
  virtual double getFarthestExtentFromOrigin() const{
    // this is a funny situation, since planes are technically infinte...
    // in order to have a sane answer, we just return the offset from the origin.
    // (multiplied by root 3, which was done in the old converter, why?)
    return sqrt(3.0) *  std::fabs(offset);
  }

protected:
  virtual iBase_EntityHandle getHandle( bool positive, iGeom_Instance& igm, double world_size){

    int igm_result;
    iBase_EntityHandle world_sphere = makeWorldSphere(igm, world_size);
    iBase_EntityHandle hemisphere;
    // note the reversal of sense in this call; mcnp and igeom define it differently.
    iGeom_sectionEnt( igm, world_sphere, 
                      normal.v[0], normal.v[1], normal.v[2], offset, !positive, &hemisphere, &igm_result);
    CHECK_IGEOM( igm_result, "Sectioning world for a plane" );
    return hemisphere;


  }

};


typedef  enum { X=0, Y=1, Z=2 } axis_t;

class GeneralQuadraticSurface : public SurfaceVolume {

protected:
  double A_,B_,C_,D_,E_,F_,G_,H_,J_,K_;
  int type;
  Vector3d translation;
  double rotation_mat[9];
  double extents[3];

public:
  GeneralQuadraticSurface(double A, double B, double C, double D, double E, double F, double G, double H, double J, double K):
    SurfaceVolume(),A_(A),B_(B),C_(C),D_(D),E_(E),F_(F),G_(G),H_(H),J_(J),K_(K)
  {}

};

class CylinderSurface : public SurfaceVolume {

protected:
  axis_t axis;
  double radius;
  Vector3d center;
  bool onaxis;

public:
  CylinderSurface( axis_t axis_p, double radius_p ):
    SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(true)
  {}

  CylinderSurface( axis_t axis_p, double radius_p, double trans1, double trans2 ):
    SurfaceVolume(), axis(axis_p), radius(radius_p), center(origin), onaxis(false)
  {
    switch(axis){
    case X: center.v[Y] += trans1; center.v[Z] += trans2; break;
    case Y: center.v[X] += trans1; center.v[Z] += trans2; break;
    case Z: center.v[X] += trans1; center.v[Y] += trans2; break;
    }
  }
  
  virtual double getFarthestExtentFromOrigin( ) const{
    return radius + center.length();
{'fraction_non_alphanumeric': 0.11854516386890487, 'fraction_numerical': 0.012070343725019984, 'mean_word_length': 3.722914307285768, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5846457', 'n_tokens_mistral': 4533, 'n_tokens_neox': 4304, 'n_words': 1205}