|
cgma
|
#include <iostream>#include <map>#include <cfloat>#include <cassert>#include "MCNPInput.hpp"#include "volumes.hpp"#include "geometry.hpp"#include "options.hpp"Go to the source code of this file.
| iBase_EntityHandle applyReverseTransform | ( | const Transform & | tx, |
| iGeom_Instance & | igm, | ||
| iBase_EntityHandle & | e | ||
| ) |
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;
}
| iBase_EntityHandle applyTransform | ( | const Transform & | t, |
| iGeom_Instance & | igm, | ||
| iBase_EntityHandle & | 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;
}
| static Transform axesImage | ( | const Vector3d & | v1, |
| const Vector3d & | v2, | ||
| const Vector3d & | v3, | ||
| const Vector3d & | translation = Vector3d() |
||
| ) | [static] |
Definition at line 463 of file volumes.cpp.
{
Vector3d x(1, 0, 0), y(0, 1, 0), z(0, 0, 1);
Vector3d a1 = v1.normalize(), a2 = v2.normalize(), a3 = v3.normalize();
if( OPT_DEBUG ) std::cout << "Axes image: " << a1 << " : " << a2 << " : " << a3 << std::endl;
double rot_matrix[9] =
{ a1.dot(x), a2.dot(x), a3.dot(x),
a1.dot(y), a2.dot(y), a3.dot(y),
a1.dot(z), a2.dot(z), a3.dot(z) };
return Transform( rot_matrix, translation );
}
| static iBase_EntityHandle embedWithinWorld | ( | bool | positive, |
| iGeom_Instance & | igm, | ||
| double | world_size, | ||
| iBase_EntityHandle | body, | ||
| bool | bound_with_world | ||
| ) | [static] |
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)
Definition at line 33 of file volumes.cpp.
{
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;
}
| static Transform imageZAxisTo | ( | const Vector3d & | v, |
| const Vector3d & | translation = Vector3d() |
||
| ) | [static] |
Definition at line 479 of file volumes.cpp.
{
// approach: find two vectors perpendicular to v, then use axesImage.
Vector3d a = v.normalize();
Vector3d x(1, 0, 0);
Vector3d b = x.cross( v.normalize() );
if( b.length() < DBL_EPSILON ){
// v is indistinguishable from the x axis
b = Vector3d(0, -1, 0);
if( OPT_DEBUG ) std::cout << "imageZAxisTo: Changing v " << std::endl;
}
Vector3d c = b.cross( v.normalize() );
return axesImage( c, b, a, translation );
}
| 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;
}
| bool sqIsEllipsoid | ( | const std::vector< double > & | args | ) |
Definition at line 785 of file volumes.cpp.
{
bool isEllipsoid = true;
if ( args.at(3)*args.at(3) + args.at(4)*args.at(4) + args.at(5)*args.at(5) != 0)
isEllipsoid = false;
if ( args.at(0)*args.at(1)*args.at(2)*args.at(6) >= 0 )
isEllipsoid = false;
return isEllipsoid;
}
VolumeCache default_volume_cache [static] |
Definition at line 800 of file volumes.cpp.
1.7.6.1