cgma
AnalyticGeometryTool Class Reference

Performs calculations on analytic geometry. More...

#include <AnalyticGeometryTool.hpp>

List of all members.

Public Member Functions

 ~AnalyticGeometryTool ()
CubitBoolean dbl_eq (double val_1, double val_2)
double get_epsilon ()
double set_epsilon (double new_epsilon)
void round_near_val (double &val)
void transform_pnt (double m[4][4], double pin[3], double pout[3])
void transform_vec (double m3[3][3], double vin[3], double vout[3])
void transform_vec (double m4[4][4], double vin[3], double vout[3])
void transform_line (double rot_mtx[4][4], double origin[3], double axis[3])
void transform_line (double rot_mtx[4][4], CubitVector &origin, CubitVector &axis)
void add_origin_to_rotation_mtx (double rot_mtx[4][4], double origin[3])
double det_mtx (double m[3][3])
 Find determinant of matrix.
CubitStatus inv_mtx_adj (double mtx[3][3], double inv_mtx[3][3])
CubitStatus inv_trans_mtx (double transf[4][4], double inv_transf[4][4])
void copy_pnt (double pnt_in[3], double pnt_out[3])
void set_pnt (double pnt[3], double x, double y, double z)
 Sets the value of pnt to x,y,z (inline function).
CubitBoolean pnt_eq (double pnt1[3], double pnt2[3])
CubitStatus mirror_pnt (double pnt[3], double pln_orig[3], double pln_norm[3], double m_pnt[3])
CubitStatus next_pnt (double str_pnt[3], double vec_dir[3], double len, double new_pnt[3])
CubitStatus unit_vec (double vin[3], double vout[3])
 Finds unit vector of input vector (vout can equal vin - converts in place).
double dot_vec (double uval[3], double vval[3])
void cross_vec (double uval[3], double vval[3], double cross[3])
void cross_vec_unit (double uval[3], double vval[3], double cross[3])
 Finds unit vector that is cross product of two vectors.
void orth_vecs (double uvect[3], double vvect[3], double wvect[3])
 Finds 2 arbitrary orthoganal vectors to the first.
double mag_vec (double vec[3])
CubitStatus get_vec (double str_pnt[3], double stp_pnt[3], double vector_out[3])
CubitStatus get_vec_unit (double str_pnt[3], double stp_pnt[3], double uv_out[3])
void mult_vecxconst (double constant, double vec[3], double vec_out[3])
 Multiples a vector by a constant (vec_out can equal vec).
void reverse_vec (double vin[3], double vout[3])
 Multiples -1.0 by a vector (vout can equal vin).
double angle_vec_vec (double v1[3], double v2[3])
double dist_pnt_pnt (double pnt1[3], double pnt2[3])
 Determines distance between two points in space.
double dist_pln_pln (double pln_1_orig[3], double pln_1_norm[3], double pln_2_orig[3], double pln_2_norm[3], AgtSide *side=NULL, AgtOrientation *orien=NULL, unsigned short *status=NULL)
CubitStatus int_ln_pln (double ln_orig[3], double ln_vec[3], double pln_orig[3], double pln_norm[3], double int_pnt[3])
int int_ln_ln (double p1[3], double v1[3], double p2[3], double v2[3], double int_pnt1[3], double int_pnt2[3])
int int_pnt_ln (double pnt[3], double ln_orig[3], double ln_vec[3], double int_pnt[3])
int int_pnt_pln (double pnt[3], double pln_orig[3], double pln_norm[3], double pln_int[3])
CubitStatus int_pln_pln (double pln_1_orig[3], double pln_1_norm[3], double pln_2_orig[3], double pln_2_norm[3], double ln_orig[3], double ln_vec[3])
CubitBoolean is_vec_par (double vec_1[3], double vec_2[3])
CubitBoolean is_vec_perp (double vec_1[3], double vec_2[3])
CubitBoolean is_vecs_same_dir (double vec_1[3], double vec_2[3])
CubitBoolean is_pnt_on_ln (double pnt[3], double ln_orig[3], double ln_vec[3])
CubitBoolean is_pnt_on_ln_seg (double pnt1[3], double end1[3], double end2[3])
CubitBoolean is_pnt_on_pln (double pnt[3], double pln_orig[3], double pln_norm[3])
CubitBoolean is_ln_on_pln (double ln_orig[3], double ln_vec[3], double pln_orig[3], double pln_norm[3])
CubitBoolean is_pln_on_pln (double pln_orig1[3], double pln_norm1[3], double pln_orig2[3], double pln_norm2[3])
 This routine checks to see if two infinite planes are coincident.
int get_num_circle_tess_pnts (double radius, double len_tolerance=1e-4)
void get_pln_orig_norm (double A, double B, double C, double D, double pln_orig[3], double pln_norm[3]=NULL)
void get_box_corners (double box_min[3], double box_max[3], double c[8][3])
double MinTriangleTriangle (Triangle3 &tri0, Triangle3 &tri1, double &s, double &t, double &u, double &v)
 Finds minimum distance between two triangles in 3D (MAGIC)
double AreaIntersection (Triangle &tri1, Triangle &tri2)
 Finds area of intersection between two triangles (MAGIC)
double Angle (Triangle3 &tri1, Triangle3 &tri2)
 Finds angle between normals of the two triangles (radians)
double MinPointTriangle (const Point3 &p, const Triangle3 &tri, double &s, double &t)
 Finds minimum distance beween a triange and point in 3D (MAGIC)
void Normal (Triangle3 &tri, double normal[3])
 Finds normal vector of given triangle.
double ProjectedOverlap (Triangle3 &tri1, Triangle3 &tri2, bool draw_overlap=false)
 Projects tri2 to the plane of tri1 and finds the overlap area.
void copy_mtx (double from[3][3], double to[3][3])
void copy_mtx (double from[4][4], double to[4][4])
void copy_mtx (double from[4][4], double to[3][3])
void copy_mtx (double from[3][3], double to[4][4], double *origin=NULL)
void create_rotation_mtx (double theta, double v[3], double mtx3x3[3][3])
void create_rotation_mtx (double theta, double v[3], double mtx4x4[4][4])
void identity_mtx (double mtx3x3[3][3])
 Simply sets the given matrix to the identity matrix.
void identity_mtx (double mtx4x4[4][4])
 Simply sets the given matrix to the identity matrix.
void mtx_to_angs (double mtx3x3[3][3], double &ax, double &ay, double &az)
void mtx_to_angs (double mtx4x4[4][4], double &ax, double &ay, double &az)
void rotate_system (double mtx[3][3], double sys_in[3][3], double sys_out[3][3])
void rotate_system (double mtx[4][4], double sys_in[4][4], double sys_out[4][4])
void mult_mtx (double a[3][3], double b[3][3], double d[3][3])
void mult_mtx (double a[4][4], double b[4][4], double d[4][4])
void vecs_to_mtx (double xvec[3], double yvec[3], double zvec[3], double matrix[3][3])
void vecs_to_mtx (double xvec[3], double yvec[3], double zvec[3], double origin[3], double matrix[4][4])
void print_mtx (double mtx[3][3])
 Prints matrix values, for debugging.
void print_mtx (double mtx[4][4])
 Prints matrix values, for debugging.
void copy_pnt (double pnt_in[3], CubitVector &cubit_vec)
 For going back and forth from CubitVectors.
void copy_pnt (CubitVector &cubit_vec, double pnt_out[3])
 For going back and forth from CubitVectors.
void setup_arc (double vec_1[3], double vec_2[3], double origin[3], double start_angle, double end_angle, double radius, AGT_3D_Arc &arc)
void setup_arc (CubitVector &vec_1, CubitVector &vec_2, CubitVector origin, double start_angle, double end_angle, double radius, AGT_3D_Arc &arc)
void get_arc_xyz (AGT_3D_Arc &arc, double param, double pnt[3])
void get_arc_xyz (AGT_3D_Arc &arc, double param, CubitVector &pnt)
CubitStatus min_pln_box_int_corners (const CubitPlane &plane, const CubitBox &box, int extension_type, double extension, CubitVector &p1, CubitVector &p2, CubitVector &p3, CubitVector &p4, CubitBoolean silent=CUBIT_FALSE)
CubitStatus min_pln_box_int_corners (CubitVector &vec1, CubitVector &vec2, CubitVector &vec3, CubitVector &box_min, CubitVector &box_max, int extension_type, double extension, CubitVector &p1, CubitVector &p2, CubitVector &p3, CubitVector &p4, CubitBoolean silent=CUBIT_FALSE)
CubitStatus min_pln_box_int_corners (double plane_norm[3], double plane_coeff, double box_min[3], double box_max[3], int extension_type, double extension, double p1[3], double p2[3], double p3[3], double p4[3], CubitBoolean silent=CUBIT_FALSE)
CubitStatus get_tight_bounding_box (DLIList< CubitVector * > &point_list, CubitVector &p1, CubitVector &p2, CubitVector &p3, CubitVector &p4, CubitVector &p5, CubitVector &p6, CubitVector &p7, CubitVector &p8)
CubitStatus get_tight_bounding_box (DLIList< CubitVector * > &point_list, CubitVector &center, CubitVector axes[3], CubitVector &extension)
CubitStatus min_cyl_box_int (double radius, CubitVector &axis, CubitVector &center, CubitBox &box, int extension_type, double extension, CubitVector &start, CubitVector &end, int num_tess_pnts=2048)
CubitStatus min_cyl_box_int (double radius, double axis[3], double center[3], double box_min[3], double box_max[3], int extension_type, double extension, double start[3], double end[3], int num_tess_pnts=2048)

Static Public Member Functions

static AnalyticGeometryToolinstance ()
static void delete_instance ()

Protected Member Functions

 AnalyticGeometryTool ()

Private Member Functions

double Area (int N, Point2 *pt, double angle)
void MinimalBoxForAngle (int N, Point2 *pt, double angle, OBBox2 &box)
OBBox2 MinimalBox2 (int N, Point2 *pt)
void MatrixToAngleAxis (double **R, double &angle, double axis[3])
void AngleAxisToMatrix (double angle, double axis[3], double R[3][3])
double Volume (int N, Point3 *pt, double angle[3])
void MinimalBoxForAngles (int N, Point3 *pt, double angle[3], OBBox3 &box)
void GetInterval (double A[3], double D[3], double &tmin, double &tmax)
void Combine (double result[3], double A[3], double t, double D[3])
double MinimizeOnInterval (int N, Point3 *pt, double A[3], double D[3])
double MinimizeOnLattice (int N, Point3 *pt, double A[3], int layers, double thickness)
void InitialGuess (int N, Point3 *pt, double angle[3])
OBBox3 MinimalBox3 (int N, Point3 *pt)
double Abs (double x)
double ACos (double x)
double ATan2 (double y, double x)
double Cos (double x)
double Sign (double x)
double Sin (double x)
double Sqrt (double x)
double UnitRandom ()
double Dot (const Point2 &p, const Point2 &q)
double Length (const Point2 &p)
CubitBoolean Unitize (Point2 &p)
double Dot (const Point3 &p, const Point3 &q)
Point3 Cross (const Point3 &p, const Point3 &q)
double Length (const Point3 &p)
CubitBoolean Unitize (Point3 &p)
double MinLineSegmentLineSegment (const Line3 &seg0, const Line3 &seg1, double &s, double &t)
int get_plane_bbox_intersections (double box_min[3], double box_max[3], double pln_orig[3], double pln_norm[3], double int_array[6][3])
double MinLineSegmentTriangle (const Line3 &seg, const Triangle3 &tri, double &r, double &s, double &t)
AgtLine EdgeToLine (Point2 *v0, Point2 *v1)
void TriangleLines (Triangle *T, AgtLine line[3])
AgtTriListSplitAndDecompose (Triangle *T, AgtLine *line)
AgtTriListIntersection (Triangle *T0, Triangle *T1)
double AreaTriangle (Triangle *T)
double Area (AgtTriList *list)
double Area (Triangle &tri1, Triangle &tri2)

Private Attributes

double agtEpsilon

Static Private Attributes

static AnalyticGeometryToolinstance_ = 0

Detailed Description

Performs calculations on analytic geometry.

Definition at line 271 of file AnalyticGeometryTool.hpp.


Constructor & Destructor Documentation

Definition at line 142 of file AnalyticGeometryTool.cpp.

{
   agtEpsilon = 1e-8;
}

Member Function Documentation

double AnalyticGeometryTool::Abs ( double  x) [inline, private]

Definition at line 965 of file AnalyticGeometryTool.hpp.

{
    return double(fabs(x));
}
double AnalyticGeometryTool::ACos ( double  x) [private]
void AnalyticGeometryTool::add_origin_to_rotation_mtx ( double  rot_mtx[4][4],
double  origin[3] 
)

Adds origin to rotation matrix, so you can rotate points about a line. Line is defined as original vector used in create_rotation_mtx and the origin.

Definition at line 412 of file AnalyticGeometryTool.cpp.

{  
   double tmp_mtx[4][4];

   // Translate to origin
   double t[4][4]; 
   memcpy(t, AGT_IDENTITY_MTX_4X4, sizeof(double)*16);
   //PRINT_INFO( "Rotation matrix, before origin: \n" );
   //print_mtx( rot_mtx );
   t[3][0]=-origin[0]; t[3][1]=-origin[1]; t[3][2]=-origin[2];
   mult_mtx( t, rot_mtx, tmp_mtx );
   
   //PRINT_INFO( "Origin times rotation: \n" );
   //print_mtx( tmp_mtx );
   
   // Translate back
   t[3][0]=origin[0]; t[3][1]=origin[1]; t[3][2]=origin[2];
   mult_mtx( tmp_mtx, t, rot_mtx );
   
   //PRINT_INFO( "Rotation x -origin: \n" );
   //print_mtx( rot_mtx );
}
double AnalyticGeometryTool::Angle ( Triangle3 tri1,
Triangle3 tri2 
)

Finds angle between normals of the two triangles (radians)

Definition at line 6204 of file AnalyticGeometryTool.cpp.

{
  double norm1[3];
  Normal( tri1, norm1 );

  double norm2[3];
  Normal( tri2, norm2 );

  return angle_vec_vec( norm1, norm2 );
}
double AnalyticGeometryTool::angle_vec_vec ( double  v1[3],
double  v2[3] 
)

Finds angle in radians between two vectors. Angle will always be between 0.0 and PI. For accuracy, the routine uses the cosine for large angles and the sine for small angles.

Definition at line 1254 of file AnalyticGeometryTool.cpp.

{   
  double denom, dot, cosang, sinang, acrsb, angle;
  double crossed_vec[3];
     
  // For accuracy, use cosine for large angles, sine for small angles 
  denom = sqrt((v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2])
              *(v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]));
  
  // Check for a zero length vector 
  if (dbl_eq(denom, 0.0))
     return (0.0);
  
  dot = dot_vec(v1, v2);
  
  cosang = dot/denom;
  
  if (1.0 - fabs(cosang) < 0.01) {
     cross_vec(v1, v2, crossed_vec);
     acrsb = mag_vec(crossed_vec);
     sinang = acrsb/denom;
     if (cosang > 0.0) 
        angle = asin(sinang);
     else 
        angle = AGT_PI - asin(sinang);
  }
  else 
     angle = acos(cosang);
  
  return (angle); 
}
void AnalyticGeometryTool::AngleAxisToMatrix ( double  angle,
double  axis[3],
double  R[3][3] 
) [private]

Definition at line 2822 of file AnalyticGeometryTool.cpp.

{
    double cs = cos(angle), sn = sin(angle);
    double length = sqrt(axis[0]*axis[0]+axis[1]*axis[1]+axis[2]*axis[2]);
    double x = axis[0]/length;
    double y = axis[1]/length;
    double z = axis[2]/length;
    double omc = 1.0-cs;
    double x2 = x*x, y2 = y*y, z2 = z*z;
    double xy = x*y, xz = x*z, yz = y*z;
    double snx = sn*x, sny = sn*y, snz = sn*z;
    
    R[0][0] = 1.0-omc*(y2+z2);
    R[0][1] = +snz+omc*xy;
    R[0][2] = -sny+omc*xz;
    R[1][0] = -snz+omc*xy;
    R[1][1] = 1.0-omc*(x2+z2);
    R[1][2] = +snx+omc*yz;
    R[2][0] = +sny+omc*xz;
    R[2][1] = -snx+omc*yz;
    R[2][2] = 1.0-omc*(x2+y2);
}
double AnalyticGeometryTool::Area ( int  N,
Point2 pt,
double  angle 
) [private]

Definition at line 2547 of file AnalyticGeometryTool.cpp.

{
    double cs = cos(angle), sn = sin(angle);

    double umin = +cs*pt[0].x+sn*pt[0].y, umax = umin;
    double vmin = -sn*pt[0].x+cs*pt[0].y, vmax = vmin;
    for (int i = 1; i < N; i++)
    {
        double u = +cs*pt[i].x+sn*pt[i].y;
        if ( u < umin )
            umin = u;
        else if ( u > umax )
            umax = u;

        double v = -sn*pt[i].x+cs*pt[i].y;
        if ( v < vmin )
            vmin = v;
        else if ( v > vmax )
            vmax = v;
    }

    double area = (umax-umin)*(vmax-vmin);
    return area;
}
double AnalyticGeometryTool::Area ( AgtTriList list) [private]

Definition at line 6145 of file AnalyticGeometryTool.cpp.

{
    double area = 0.0;

    while ( list )
    {
        area += AreaTriangle(list->tri);
        list = list->next;
    }

    return area;
}
double AnalyticGeometryTool::Area ( Triangle tri1,
Triangle tri2 
) [private]
double AnalyticGeometryTool::AreaIntersection ( Triangle tri1,
Triangle tri2 
)

Finds area of intersection between two triangles (MAGIC)

Definition at line 6160 of file AnalyticGeometryTool.cpp.

{
  AgtTriList* list = Intersection(&tri1,&tri2);
  double area = Area(list);
  while ( list )
  {
    AgtTriList* save = list;
    list = list->next;
    delete save->tri;
    delete save;
  }
  delete list;
  return area;
}
double AnalyticGeometryTool::AreaTriangle ( Triangle T) [private]

Definition at line 6136 of file AnalyticGeometryTool.cpp.

{
    Point2 e1 = { T->v[1].x - T->v[0].x, T->v[1].y - T->v[0].y };
    Point2 e2 = { T->v[2].x - T->v[0].x, T->v[2].y - T->v[0].y };

    return fabs(0.5*(e1.x*e2.y-e1.y*e2.x));
}
double AnalyticGeometryTool::ATan2 ( double  y,
double  x 
) [inline, private]

Definition at line 970 of file AnalyticGeometryTool.hpp.

{
    return double(atan2(y,x));
}
void AnalyticGeometryTool::Combine ( double  result[3],
double  A[3],
double  t,
double  D[3] 
) [private]

Definition at line 3009 of file AnalyticGeometryTool.cpp.

{
    for (int i = 0; i < 3; i++)
        result[i] = A[i]+t*D[i];
}
void AnalyticGeometryTool::copy_mtx ( double  from[3][3],
double  to[3][3] 
)

This routine simply copies one matrix to another. If a NULL is passed in for the from matrix, then the identity matrix is copied into the out matrix. Last function allows you to specify 1st 3 doubles of row 4 (origin) of the 4x4 matrix.

Definition at line 271 of file AnalyticGeometryTool.cpp.

{   
   // Determine if identity matrix needed 
   if (!from) // copy in the identity matrix 
      memcpy(to, AGT_IDENTITY_MTX_3X3, sizeof(double)*9); 
   else // Copy from to to 
      memcpy(to, from, sizeof(double)*9);   
}
void AnalyticGeometryTool::copy_mtx ( double  from[4][4],
double  to[4][4] 
)

This routine simply copies one matrix to another. If a NULL is passed in for the from matrix, then the identity matrix is copied into the out matrix. Last function allows you to specify 1st 3 doubles of row 4 (origin) of the 4x4 matrix.

Definition at line 280 of file AnalyticGeometryTool.cpp.

{   
   // determine if identity matrix needed 
   if (!from) // copy in the identity matrix 
      memcpy(to, AGT_IDENTITY_MTX_4X4, sizeof(double)*16);  
   else // copy from to to 
      memcpy(to, from, sizeof(double)*16);   
}
void AnalyticGeometryTool::copy_mtx ( double  from[4][4],
double  to[3][3] 
)

This routine simply copies one matrix to another. If a NULL is passed in for the from matrix, then the identity matrix is copied into the out matrix. Last function allows you to specify 1st 3 doubles of row 4 (origin) of the 4x4 matrix.

Definition at line 289 of file AnalyticGeometryTool.cpp.

{   
   size_t dbl3;
   
   dbl3 = sizeof(double) * 3;
   
   // Determine if identity matrix needed 
   if (!from) // Copy in the identity matrix 
      memcpy(to, AGT_IDENTITY_MTX_3X3, sizeof(double)*9); 
   else { // Copy each upper left element of from to to 
      memcpy(to[0], from[0], dbl3);
      memcpy(to[1], from[1], dbl3);
      memcpy(to[2], from[2], dbl3);
   }   
}
void AnalyticGeometryTool::copy_mtx ( double  from[3][3],
double  to[4][4],
double *  origin = NULL 
)

This routine simply copies one matrix to another. If a NULL is passed in for the from matrix, then the identity matrix is copied into the out matrix. Last function allows you to specify 1st 3 doubles of row 4 (origin) of the 4x4 matrix.

Definition at line 305 of file AnalyticGeometryTool.cpp.

{   
   size_t dbl3;
   
   dbl3 = sizeof(double) * 3;
   
   // Determine if identity matrix needed 
   if (!from) { // Copy in the identity matrix 
   
      memcpy(to, AGT_IDENTITY_MTX_4X4, sizeof(double)*16);
      
      if (origin)
         memcpy(to[3], origin, dbl3);
   }   
         
   else { // Copy each upper element of from to to 
      
      memcpy(to[0], from[0], dbl3);
      memcpy(to[1], from[1], dbl3);
      memcpy(to[2], from[2], dbl3);
      
      to[0][3] = 0.0;
      to[1][3] = 0.0;
      to[2][3] = 0.0;
      to[3][3] = 1.0;
      
      if (origin)
      {
         memcpy(to[3], origin, dbl3);
//         to[3][0] = origin[0];
//         to[3][1] = origin[1];
//         to[3][2] = origin[2];
      }
      else
         memcpy(to[3], AGT_IDENTITY_MTX_4X4[3], sizeof(double)*3);
   } 
}
void AnalyticGeometryTool::copy_pnt ( double  pnt_in[3],
double  pnt_out[3] 
)

Copy one double[3] point to another double[3] point (uses memcpy). If first point in NULL then second point is set to 0,0,0.

Definition at line 998 of file AnalyticGeometryTool.cpp.

{   
   if (pnt_in == pnt_out)
      return;
   
   if (pnt_out == NULL)
      return;
   
   if (pnt_in == NULL) {
      pnt_out[0] = 0.0;
      pnt_out[1] = 0.0;
      pnt_out[2] = 0.0;
      return;
   }
   
   // Simply copy first point into second point    
   memcpy(pnt_out, pnt_in, sizeof(double)*3);
}
void AnalyticGeometryTool::copy_pnt ( double  pnt_in[3],
CubitVector cubit_vec 
)

For going back and forth from CubitVectors.

Definition at line 1017 of file AnalyticGeometryTool.cpp.

{
   cubit_vec.set( pnt_in );
}
void AnalyticGeometryTool::copy_pnt ( CubitVector cubit_vec,
double  pnt_out[3] 
)

For going back and forth from CubitVectors.

Definition at line 1022 of file AnalyticGeometryTool.cpp.

{
   pnt_out[0] = cubit_vec.x();
   pnt_out[1] = cubit_vec.y();
   pnt_out[2] = cubit_vec.z();
}
double AnalyticGeometryTool::Cos ( double  x) [private]
void AnalyticGeometryTool::create_rotation_mtx ( double  theta,
double  v[3],
double  mtx3x3[3][3] 
)

This routine determines the tranformation matrix given the theta and the vector to cross through.

Definition at line 344 of file AnalyticGeometryTool.cpp.

{        
   double coeff1;
   double coeff2;
   double coeff3;
   double v_unit[3];
   
   if (!mtx3x3)
      return;
   
   coeff1 = cos(theta);
   coeff2 = (1.0l - coeff1);
   coeff3 = sin(theta);
   
   unit_vec(v, v_unit);
   
   mtx3x3[0][0] = coeff1 + coeff2 * (v_unit[0] * v_unit[0]);
   mtx3x3[0][1] = coeff2 * v_unit[1] * v_unit[0] + coeff3 * v_unit[2];
   mtx3x3[0][2] = coeff2 * v_unit[2] * v_unit[0] - coeff3 * v_unit[1];
   
   mtx3x3[1][0] = coeff2 * v_unit[1] * v_unit[0] - coeff3 * v_unit[2];
   mtx3x3[1][1] = coeff1 + coeff2 * (v_unit[1] * v_unit[1]);
   mtx3x3[1][2] = coeff2 * v_unit[1] * v_unit[2] + coeff3 * v_unit[0];
   
   mtx3x3[2][0] = coeff2 * v_unit[2] * v_unit[0] + coeff3 * v_unit[1];
   mtx3x3[2][1] = coeff2 * v_unit[2] * v_unit[1] - coeff3 * v_unit[0];
   mtx3x3[2][2] = coeff1 + coeff2 * (v_unit[2] * v_unit[2]);
}
void AnalyticGeometryTool::create_rotation_mtx ( double  theta,
double  v[3],
double  mtx4x4[4][4] 
)

This routine determines the tranformation matrix given the theta and the vector to cross through.

Definition at line 374 of file AnalyticGeometryTool.cpp.

{        
   double coeff1;
   double coeff2;
   double coeff3;
   double v_unit[3];
   
   if (!mtx4x4)
      return;
   
   coeff1 = cos(theta);
   coeff2 = (1.0l - coeff1);
   coeff3 = sin(theta);

   unit_vec(v, v_unit);
   
   mtx4x4[0][0] = coeff1 + coeff2 * (v_unit[0] * v_unit[0]);
   mtx4x4[0][1] = coeff2 * v_unit[1] * v_unit[0] + coeff3 * v_unit[2];
   mtx4x4[0][2] = coeff2 * v_unit[2] * v_unit[0] - coeff3 * v_unit[1];
   
   mtx4x4[1][0] = coeff2 * v_unit[1] * v_unit[0] - coeff3 * v_unit[2];
   mtx4x4[1][1] = coeff1 + coeff2 * (v_unit[1] * v_unit[1]);
   mtx4x4[1][2] = coeff2 * v_unit[1] * v_unit[2] + coeff3 * v_unit[0];
   
   mtx4x4[2][0] = coeff2 * v_unit[2] * v_unit[0] + coeff3 * v_unit[1];
   mtx4x4[2][1] = coeff2 * v_unit[2] * v_unit[1] - coeff3 * v_unit[0];
   mtx4x4[2][2] = coeff1 + coeff2 * (v_unit[2] * v_unit[2]);
   
   mtx4x4[0][3] = 0.0;
   mtx4x4[1][3] = 0.0;
   mtx4x4[2][3] = 0.0;
   mtx4x4[3][3] = 1.0;
   mtx4x4[3][0] = 0.0;
   mtx4x4[3][1] = 0.0;
   mtx4x4[3][2] = 0.0;  
}
Point3 AnalyticGeometryTool::Cross ( const Point3 p,
const Point3 q 
) [private]
void AnalyticGeometryTool::cross_vec ( double  uval[3],
double  vval[3],
double  cross[3] 
)

Finds cross product of two vectors: cross[0] = u[1] * v[2] - u[2] * v[1]; cross[1] = u[2] * v[0] - u[0] * v[2]; cross[2] = u[0] * v[1] - u[1] * v[0];

Definition at line 1128 of file AnalyticGeometryTool.cpp.

{      
   // determine cross product of the two vectors 
   
   cross[0] = uval[1] * vval[2] - uval[2] * vval[1];
   cross[1] = uval[2] * vval[0] - uval[0] * vval[2];
   cross[2] = uval[0] * vval[1] - uval[1] * vval[0];   
}
void AnalyticGeometryTool::cross_vec_unit ( double  uval[3],
double  vval[3],
double  cross[3] 
)

Finds unit vector that is cross product of two vectors.

Definition at line 1138 of file AnalyticGeometryTool.cpp.

{      
   // determine cross product of the two vectors 
   cross_vec(uval,vval,cross);
   
   // convert to unit vector 
   unit_vec(cross,cross);
}
CubitBoolean AnalyticGeometryTool::dbl_eq ( double  val_1,
double  val_2 
) [inline]

Check to see if double numbers are equal within epsilon. Values are equal if fabs(val_1 - val_2) < epsilon. Epsilon is determined by next two functions.

Definition at line 1002 of file AnalyticGeometryTool.hpp.

{ 
   CubitBoolean result;
   if (fabs(val_1 - val_2) < agtEpsilon)
      result = CUBIT_TRUE;
   else
      result = CUBIT_FALSE;
   return result;
}
static void AnalyticGeometryTool::delete_instance ( ) [inline, static]

Definition at line 278 of file AnalyticGeometryTool.hpp.

   {
     if( instance_ )
       delete instance_;
     instance_ = NULL;
   };
double AnalyticGeometryTool::det_mtx ( double  m[3][3])

Find determinant of matrix.

Definition at line 693 of file AnalyticGeometryTool.cpp.

{
   double determinant;
   
   if (!m)
      return (0.0);
   
   determinant = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
               + m[0][1]*(m[1][2]*m[2][0]-m[1][0]*m[2][2])
               + m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
   
   return (determinant);
}  
double AnalyticGeometryTool::dist_pln_pln ( double  pln_1_orig[3],
double  pln_1_norm[3],
double  pln_2_orig[3],
double  pln_2_norm[3],
AgtSide side = NULL,
AgtOrientation orien = NULL,
unsigned short *  status = NULL 
)

This routine determines the shortest distance between two planes. This is the perpendicular distance between the two planes. Note that if the two planes are not parallel, the returned distance is zero. The routine can also be used to determine which side of the first plane the second plane is on, and the orientation of the two planes to each other.

Definition at line 1301 of file AnalyticGeometryTool.cpp.

{
   double distance;
   double int_pnt[3];
   double vec[3];
   
   // Check to see if planes are parallel
   if (is_vec_par(pln_1_norm, pln_2_norm)) {    
      
      // Set successful status
      if (status)
         *status = 1;
      
      // Calculate perpendicular line plane intersection on plane_2 from 
      // pln_1_origin 
      int_ln_pln(pln_1_orig, pln_1_norm, pln_2_orig, pln_2_norm, 
                 int_pnt);
      
      // Find distance between pln_1_origin and this intersection pnt
      distance = dist_pnt_pnt(pln_1_orig, int_pnt);
      
      // Get side if required
      if (side) {
         
         if (dbl_eq(distance, 0.0))
            
            *side = AGT_ON_PLANE;
         
         else {
            
            // Get vector to intersection point
            get_vec(pln_1_orig, int_pnt, vec);
            
            // Compare angles
            if (dbl_eq(angle_vec_vec(vec, pln_1_norm), 0.0))
               *side = AGT_POS_SIDE;
            else
               *side = AGT_NEG_SIDE;
            
         }
      }
      
      // Get orientation if required
      if (orien) {
         
         // Compare surface normals
         if (dbl_eq(angle_vec_vec(pln_1_norm, pln_2_norm), 0.0))
            *orien = AGT_SAME_DIR;
         else
            *orien = AGT_OPP_DIR;
      }      
      
   }
   
   else {
      
      if (status)
         *status = 0;
      
      if (side)
         *side = AGT_CROSS;
      
      distance = 0.0;
      
   }
   
   return (distance);
}
double AnalyticGeometryTool::dist_pnt_pnt ( double  pnt1[3],
double  pnt2[3] 
)

Determines distance between two points in space.

Definition at line 1289 of file AnalyticGeometryTool.cpp.

{  
  double x = pnt2[0] - pnt1[0];  // difference in the x direction
  double y = pnt2[1] - pnt1[1];  // difference in the y direction
  double z = pnt2[2] - pnt1[2];  // difference in the z direction
  
  // return the distance   
  return(sqrt(x*x + y*y + z*z));
}
double AnalyticGeometryTool::Dot ( const Point2 p,
const Point2 q 
) [inline, private]

Definition at line 81 of file AnalyticGeometryTool.cpp.

{
    return double(p.x*q.x + p.y*q.y);
}
double AnalyticGeometryTool::Dot ( const Point3 p,
const Point3 q 
) [inline, private]

Definition at line 989 of file AnalyticGeometryTool.hpp.

{
    return double(p.x*q.x + p.y*q.y + p.z*q.z);
}
double AnalyticGeometryTool::dot_vec ( double  uval[3],
double  vval[3] 
)

Dots two vectors. Property of dot product is: angle between vectors acute if dot product > 0 angle between vectors obtuse if dot product < 0 angle between vectors 90 deg if dot product = 0

Definition at line 1117 of file AnalyticGeometryTool.cpp.

{
   double dot_val;
   
   // perform dot calculation = v[0]*u[0] + v[1]*u[1] + v[1]*u[1] 
        
   dot_val = uval[0]*vval[0] + uval[1]*vval[1] + uval[2]*vval[2];
      
   return(dot_val);   
}
AgtLine AnalyticGeometryTool::EdgeToLine ( Point2 v0,
Point2 v1 
) [private]

Definition at line 5942 of file AnalyticGeometryTool.cpp.

{
    // assert:  v0 and v1 are distinct

    Point2 edge = { v1->x - v0->x, v1->y - v0->y };

    AgtLine line;
    line.N.x = edge.y;
    line.N.y = -edge.x;
    line.c = line.N.x*v0->x + line.N.y*v0->y;

    return line;
}
void AnalyticGeometryTool::get_arc_xyz ( AGT_3D_Arc arc,
double  param,
double  pnt[3] 
)

Given a parameter value from 0 to 1 on the arc, return the xyz location. Arc is assumed to be parameterized from 0 to 1.

Definition at line 1761 of file AnalyticGeometryTool.cpp.

{
   double Tp;
   
   // Un-normalized parameter
   Tp = arc.StartAngle * ( 1.0 - param ) + arc.EndAngle * param;
   
   // Solve for XYZ
   pnt[0] = arc.Radius * ( cos( Tp ) * arc.Vec1[0] +
                           sin( Tp ) * arc.Vec2[0] ) + 
                           arc.Origin[0];
   
   pnt[1] = arc.Radius * ( cos( Tp ) * arc.Vec1[1] +
                           sin( Tp ) * arc.Vec2[1] ) + 
                           arc.Origin[1];
   
   pnt[2] = arc.Radius * ( cos( Tp ) * arc.Vec1[2] +
                           sin( Tp ) * arc.Vec2[2] ) + 
                           arc.Origin[2];
}
void AnalyticGeometryTool::get_arc_xyz ( AGT_3D_Arc arc,
double  param,
CubitVector pnt 
)

Given a parameter value from 0 to 1 on the arc, return the xyz location. Arc is assumed to be parameterized from 0 to 1.

Definition at line 1782 of file AnalyticGeometryTool.cpp.

{
   double Tp;
   
   // Un-normalized parameter
   Tp = arc.StartAngle * ( 1.0 - param ) + arc.EndAngle * param;
   
   // Solve for XYZ
   pnt.x( arc.Radius * ( cos( Tp ) * arc.Vec1[0] +
                         sin( Tp ) * arc.Vec2[0] ) + 
                         arc.Origin[0] );
   
   pnt.y( arc.Radius * ( cos( Tp ) * arc.Vec1[1] +
                         sin( Tp ) * arc.Vec2[1] ) + 
                         arc.Origin[1] );
   
   pnt.z( arc.Radius * ( cos( Tp ) * arc.Vec1[2] +
                         sin( Tp ) * arc.Vec2[2] ) + 
                         arc.Origin[2] );
}
void AnalyticGeometryTool::get_box_corners ( double  box_min[3],
double  box_max[3],
double  c[8][3] 
)

Find 8 corners of a box given minimum and maximum points. Box is defined starting from left-bottom-front clockwise (at front of box), then to the rear in same order.

Definition at line 1885 of file AnalyticGeometryTool.cpp.

{
   // Left-Bottom-Front       // Left-Top-Front
   c[0][0] = box_min[0];      c[1][0] = box_min[0];
   c[0][1] = box_min[1];      c[1][1] = box_max[1];
   c[0][2] = box_max[2];      c[1][2] = box_max[2];

   // Right-Top-Front         // Right-Bottom-Front
   c[2][0] = box_max[0];      c[3][0] = box_max[0];
   c[2][1] = box_max[1];      c[3][1] = box_min[1];
   c[2][2] = box_max[2];      c[3][2] = box_max[2];

   // Left-Bottom-Back        // Left-Top-Back
   c[4][0] = box_min[0];      c[5][0] = box_min[0];
   c[4][1] = box_min[1];      c[5][1] = box_max[1];
   c[4][2] = box_min[2];      c[5][2] = box_min[2];

   // Right-Top-Back          // Right-Bottom-Back
   c[6][0] = box_max[0];      c[7][0] = box_max[0];
   c[6][1] = box_max[1];      c[7][1] = box_min[1];
   c[6][2] = box_min[2];      c[7][2] = box_min[2];
   
}

Get epsilon used for double number comparisons. Default value is 1e-6.

int AnalyticGeometryTool::get_num_circle_tess_pnts ( double  radius,
double  len_tolerance = 1e-4 
)

Get number of tessellation points on the circle required to approximate the length of the circle within len_tolerance. Can be used to find the number of tessellations points to display a circle - smaller circles will have fewer tessellation points, larger circles will have more tessellation points with the same tolerance. Minimum number of tessellations points returned is 8.

Definition at line 1804 of file AnalyticGeometryTool.cpp.

{
  double cmin, cmax;

  double c = 2*CUBIT_PI*radius; // Circumference

  // Find the number of points required for the given accuracy.  Use
  // a bisection method.
  int nmin = 8, nmax = 100;
  cmin = 2.0*nmin*radius*sin(CUBIT_PI/nmin); // Circumference of circle using segments
  cmax = 2.0*nmax*radius*sin(CUBIT_PI/nmax);

  if( dbl_eq( cmin, c ) )
    return nmin;

  double old_epsilon = set_epsilon( len_tol );

  // Find an n that is more than accurate enough
  while( !dbl_eq( cmax, c ) )
  {
    nmin = nmax;
    nmax = nmin * 10;
    cmin = 2.0*nmin*radius*sin(CUBIT_PI/nmin);
    cmax = 2.0*nmax*radius*sin(CUBIT_PI/nmax);
  }
  
  // Biscect until the minimum number of segments satisfying
  // the tolerance is found.
  int n;
  while( 1 )
  {
    n = (nmin + nmax)/2;
    double cn = 2.0*n*radius*sin(CUBIT_PI/n);
    if( dbl_eq( cn, c ) )
    {
      // Go lower
      nmax = n;
    }
    else
    {
      // Go higher
      nmin = n;
    }
    if( nmax-nmin < 2 )
      break;
  }
  set_epsilon( old_epsilon );

  return nmax;
}
int AnalyticGeometryTool::get_plane_bbox_intersections ( double  box_min[3],
double  box_max[3],
double  pln_orig[3],
double  pln_norm[3],
double  int_array[6][3] 
) [private]

Definition at line 2189 of file AnalyticGeometryTool.cpp.

{
   // Fill in an array with all 8 box corners
   double corner[8][3];
   get_box_corners( box_min, box_max, corner );
 
   // Get 12 edges of the box
   double ln_start[12][3], ln_end[12][3];
   copy_pnt( corner[0], ln_start[0] );  copy_pnt( corner[1], ln_end[0] );
   copy_pnt( corner[1], ln_start[1] );  copy_pnt( corner[2], ln_end[1] );
   copy_pnt( corner[2], ln_start[2] );  copy_pnt( corner[3], ln_end[2] );
   copy_pnt( corner[3], ln_start[3] );  copy_pnt( corner[0], ln_end[3] ); 
   copy_pnt( corner[4], ln_start[4] );  copy_pnt( corner[5], ln_end[4] );
   copy_pnt( corner[5], ln_start[5] );  copy_pnt( corner[6], ln_end[5] );
   copy_pnt( corner[6], ln_start[6] );  copy_pnt( corner[7], ln_end[6] );
   copy_pnt( corner[7], ln_start[7] );  copy_pnt( corner[4], ln_end[7] );
   copy_pnt( corner[0], ln_start[8] );  copy_pnt( corner[4], ln_end[8] );
   copy_pnt( corner[1], ln_start[9] );  copy_pnt( corner[5], ln_end[9] );
   copy_pnt( corner[2], ln_start[10] ); copy_pnt( corner[6], ln_end[10] );
   copy_pnt( corner[3], ln_start[11] ); copy_pnt( corner[7], ln_end[11] );
   
   double ln_vec[3];
   double int_pnt[3];
   int num_int = 0;
   int i, j, found;
   for( i=0; i<12; i++ )
   {
      get_vec_unit( ln_start[i], ln_end[i], ln_vec );
      if( int_ln_pln( ln_start[i], ln_vec, pln_orig, pln_norm, int_pnt ) )
      {
         // Only add if on the bounded line segment
         if( is_pnt_on_ln_seg( int_pnt, ln_start[i], ln_end[i] ) )
         {
            // Only add if unique
            found = 0;
            for( j=0; j<num_int; j++ )
            {
               if( pnt_eq( int_pnt, int_array[j] ) )
               {
                  found = 1;
                  break;
               }
            }
            if( !found )
            {
               copy_pnt( int_pnt, int_array[num_int] );
               num_int++;
            }
         }
      }
   }
  return num_int;
}
void AnalyticGeometryTool::get_pln_orig_norm ( double  A,
double  B,
double  C,
double  D,
double  pln_orig[3],
double  pln_norm[3] = NULL 
)

Finds an origin-normal format plane from reduced form A*x + B*y + C*z + D = 0

Definition at line 1859 of file AnalyticGeometryTool.cpp.

{
   double x = 0.0, y = 0.0, z = 0.0;

   // Try to have origin aligned with one of the principal axes
   if( !dbl_eq( C, 0.0 ) )
      z = -D/C;
   else if (!dbl_eq( A, 0.0 ) )
      x = -D/A;
   else if (!dbl_eq( B, 0.0 ) )
      y = -D/B;

   pln_orig[0] = x;
   pln_orig[1] = y;
   pln_orig[2] = z;

   if( pln_norm )
   {
      pln_norm[0] = A;
      pln_norm[1] = B;
      pln_norm[2] = C;
   }
}

Finds the minimum size bounding box that fits around the points. Box will not necessarily be oriented in xyz directions.

Definition at line 2291 of file AnalyticGeometryTool.cpp.

{
   int num_pnts = point_list.size();
   if( num_pnts == 0 )
      return CUBIT_FAILURE;
   Point3 *pnt_arr = new Point3[num_pnts];

   int i;
   point_list.reset();
   CubitVector *vec;
   for( i=0; i<num_pnts; i++ )
   {
      vec = point_list.get_and_step();

      pnt_arr[i].x = vec->x();
      pnt_arr[i].y = vec->y();
      pnt_arr[i].z = vec->z();
   }

   OBBox3 minimal = MinimalBox3 (point_list.size(), pnt_arr);

   //PRINT_INFO( "MinBox center = %lf, %lf, %lf\n", minimal.center.x, minimal.center.y, minimal.center.z );
   //PRINT_INFO( "MinBox axis 1 = %lf, %lf, %lf\n", minimal.axis[0].x, minimal.axis[0].y, minimal.axis[0].z );
   //PRINT_INFO( "MinBox axis 2 = %lf, %lf, %lf\n", minimal.axis[1].x, minimal.axis[1].y, minimal.axis[1].z );
   //PRINT_INFO( "MinBox axis 3 = %lf, %lf, %lf\n", minimal.axis[2].x, minimal.axis[2].y, minimal.axis[2].z );
   //PRINT_INFO( "MinBox extent = %lf, %lf, %lf\n", minimal.extent[0], minimal.extent[1], minimal.extent[2] );

   CubitVector center(minimal.center.x, minimal.center.y, minimal.center.z);
   CubitVector x(minimal.axis[0].x, minimal.axis[0].y, minimal.axis[0].z);
   CubitVector y(minimal.axis[1].x, minimal.axis[1].y, minimal.axis[1].z);
   CubitVector z(minimal.axis[2].x, minimal.axis[2].y, minimal.axis[2].z);
   CubitVector extent(minimal.extent[0], minimal.extent[1], minimal.extent[2]);

   center.next_point( -x, extent.x(), p1 ); p1.next_point( -y, extent.y(), p1 );
   p1.next_point( z, extent.z(), p1 );

   center.next_point( -x, extent.x(), p2 ); p2.next_point( y, extent.y(), p2 );
   p2.next_point( z, extent.z(), p2 );

   center.next_point( x, extent.x(), p3 ); p3.next_point( y, extent.y(), p3 );
   p3.next_point( z, extent.z(), p3 );

   center.next_point( x, extent.x(), p4 ); p4.next_point( -y, extent.y(), p4 );
   p4.next_point( z, extent.z(), p4 );

   center.next_point( -x, extent.x(), p5 ); p5.next_point( -y, extent.y(), p5 );
   p5.next_point( -z, extent.z(), p5 );

   center.next_point( -x, extent.x(), p6 ); p6.next_point( y, extent.y(), p6 );
   p6.next_point( -z, extent.z(), p6 );

   center.next_point( x, extent.x(), p7 ); p7.next_point( y, extent.y(), p7 );
   p7.next_point( -z, extent.z(), p7 );

   center.next_point( x, extent.x(), p8 ); p8.next_point( -y, extent.y(), p8 );
   p8.next_point( -z, extent.z(), p8 );

   delete [] pnt_arr;

   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::get_tight_bounding_box ( DLIList< CubitVector * > &  point_list,
CubitVector center,
CubitVector  axes[3],
CubitVector extension 
)

Finds the minimum size bounding box that fits around the points. Box will not necessarily be oriented in xyz directions.

Definition at line 2248 of file AnalyticGeometryTool.cpp.

{
   int num_pnts = point_list.size();
   if( num_pnts == 0 )
      return CUBIT_FAILURE;
   Point3 *pnt_arr = new Point3[num_pnts];

   int i;
   point_list.reset();
   CubitVector *vec;
   for( i=0; i<num_pnts; i++ )
   {
      vec = point_list.get_and_step();
      //GfxPreview::draw_point( vec, 3 );

      pnt_arr[i].x = vec->x();
      pnt_arr[i].y = vec->y();
      pnt_arr[i].z = vec->z();
   }

   OBBox3 minimal = MinimalBox3 (point_list.size(), pnt_arr);

   //PRINT_INFO( "MinBox center = %lf, %lf, %lf\n", minimal.center.x, minimal.center.y, minimal.center.z );
   //PRINT_INFO( "MinBox axis 1 = %lf, %lf, %lf\n", minimal.axis[0].x, minimal.axis[0].y, minimal.axis[0].z );
   //PRINT_INFO( "MinBox axis 2 = %lf, %lf, %lf\n", minimal.axis[1].x, minimal.axis[1].y, minimal.axis[1].z );
   //PRINT_INFO( "MinBox axis 3 = %lf, %lf, %lf\n", minimal.axis[2].x, minimal.axis[2].y, minimal.axis[2].z );
   //PRINT_INFO( "MinBox extent = %lf, %lf, %lf\n", minimal.extent[0], minimal.extent[1], minimal.extent[2] );

   center.set(minimal.center.x, minimal.center.y, minimal.center.z);
   axes[0].set(minimal.axis[0].x, minimal.axis[0].y, minimal.axis[0].z);
   axes[1].set(minimal.axis[1].x, minimal.axis[1].y, minimal.axis[1].z);
   axes[2].set(minimal.axis[2].x, minimal.axis[2].y, minimal.axis[2].z);
   extension.set(minimal.extent[0], minimal.extent[1], minimal.extent[2]);

   delete [] pnt_arr;

   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::get_vec ( double  str_pnt[3],
double  stp_pnt[3],
double  vector_out[3] 
)

Finds a vector from str_pnt to stp_pnt. Returns CUBIT_FAILURE if points are coincident.

Definition at line 1201 of file AnalyticGeometryTool.cpp.

{  
   // Make sure we can create a vector 
   if (pnt_eq(str_pnt, stp_pnt)) {
      vector_out[0] = 0.0;
      vector_out[1] = 0.0;
      vector_out[2] = 0.0;
      return CUBIT_FAILURE;
   }
   
   // determine vector by subtracting starting point from stopping point 
   
  vector_out[0] = stp_pnt[0] - str_pnt[0];
  vector_out[1] = stp_pnt[1] - str_pnt[1];
  vector_out[2] = stp_pnt[2] - str_pnt[2];  
  
  return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::get_vec_unit ( double  str_pnt[3],
double  stp_pnt[3],
double  uv_out[3] 
)

Finds a unit vector from str_pnt to stp_pnt. Returns CUBIT_FAILURE if points are coincident.

Definition at line 1222 of file AnalyticGeometryTool.cpp.

{
  // determine vector between points 
   if (!get_vec(str_pnt,stp_pnt,uv_out))
      return CUBIT_FAILURE;
  
  // unitize vector 
  if (!unit_vec(uv_out,uv_out))
     return CUBIT_FAILURE;
     
  return CUBIT_SUCCESS;  
}
void AnalyticGeometryTool::GetInterval ( double  A[3],
double  D[3],
double &  tmin,
double &  tmax 
) [private]

Definition at line 2956 of file AnalyticGeometryTool.cpp.

{
    //static const double angle_min[3] = { -AGT_PI, 0.0, 0.0 };
    //static const double angle_max[3] = {  AGT_PI,  AGT_PI,  AGT_PI };
    //The pgCC compiler running on solars and cross-compiling for
    //janus has a bug such that it dies if the initialization
    //of angle_min is mixed symbolic and literal constants, so
    //define a symbolid constant for 0.0.
    // -- J.Kraftcheck, 06/05/2001
    static const double ZERO = 0.0;
    static const double angle_min[3] = { -AGT_PI, ZERO, ZERO };
    static const double angle_max[3] = {  AGT_PI, AGT_PI, AGT_PI };

    tmin = -DBL_MAX;
    tmax = +DBL_MAX;

    for (int i = 0; i < 3; i++)
    {
        const double epsilon = 1e-08;
        double b0 = angle_min[i]-A[i];
        double b1 = angle_max[i]-A[i];

        double inv, tmp;
        if ( D[i] > epsilon )
        {
            inv = 1.0/D[i];
            tmp = inv*b0;
            if ( tmp > tmin )
                tmin = tmp;
            tmp = inv*b1;
            if ( tmp < tmax )
                tmax = tmp;
        }
        else if ( D[i] < -epsilon )
        {
            inv = 1.0/D[i];
            tmp = inv*b0;
            if ( tmp < tmax )
                tmax = tmp;
            tmp = inv*b1;
            if ( tmp > tmin )
                tmin = tmp;
        }
    }

    if( tmin == -DBL_MAX || tmax == DBL_MAX ) // Added by SRS 10-6-2000
    { 
       //PRINT_WARNING( "tmin/tmax not set\n" );
       tmin = -1.0;
       tmax = 1.0;
    }
}
void AnalyticGeometryTool::identity_mtx ( double  mtx3x3[3][3])

Simply sets the given matrix to the identity matrix.

Definition at line 436 of file AnalyticGeometryTool.cpp.

{ 
   memcpy(mtx3x3, AGT_IDENTITY_MTX_3X3, sizeof(double)*9);
}
void AnalyticGeometryTool::identity_mtx ( double  mtx4x4[4][4])

Simply sets the given matrix to the identity matrix.

Definition at line 441 of file AnalyticGeometryTool.cpp.

{ 
   memcpy(mtx4x4, AGT_IDENTITY_MTX_4X4, sizeof(double)*16);
}
void AnalyticGeometryTool::InitialGuess ( int  N,
Point3 pt,
double  angle[3] 
) [private]

Definition at line 3171 of file AnalyticGeometryTool.cpp.

{
    int i;

    // compute mean of points
    double xsum = 0.0f, ysum = 0.0f, zsum = 0.0f;;
    for (i = 0; i < N; i++)
    {
        xsum += pt[i].x;
        ysum += pt[i].y;
        zsum += pt[i].z;
    }
    double xmean = xsum/N;
    double ymean = ysum/N;
    double zmean = zsum/N;

    // compute covariances of points
    double xxsum = 0.0f, xysum = 0.0f, xzsum = 0.0f;
    double yysum = 0.0f, yzsum = 0.0f, zzsum = 0.0f;
    for (i = 0; i < N; i++)
    {
        double dx = pt[i].x - xmean;
        double dy = pt[i].y - ymean;
        double dz = pt[i].z - zmean;
        xxsum += dx*dx;
        xysum += dx*dy;
        xzsum += dx*dz;
        yysum += dy*dy;
        yzsum += dy*dz;
        zzsum += dz*dz;
    }
    double xxcov = xxsum/N;
    double xycov = xysum/N;
    double xzcov = xzsum/N;
    double yycov = yysum/N;
    double yzcov = yzsum/N;
    double zzcov = zzsum/N;

    // compute eigenvectors for covariance matrix
    mgcEigenD eig(3);
    eig.Matrix(0,0) = xxcov;
    eig.Matrix(0,1) = xycov;
    eig.Matrix(0,2) = xzcov;
    eig.Matrix(1,0) = xycov;
    eig.Matrix(1,1) = yycov;
    eig.Matrix(1,2) = yzcov;
    eig.Matrix(2,0) = xzcov;
    eig.Matrix(2,1) = yzcov;
    eig.Matrix(2,2) = zzcov;
    eig.EigenStuff3();

    // Use eigenvectors as the box axes.  Eigenmatrix must not have a
    // reflection component, thus the check for negative determinant.
    const double epsilon = 1e-06;
    double** R = (double**)eig.Eigenvector();
    double det =
        +R[0][0]*R[1][1]*R[2][2]
        +R[0][1]*R[1][2]*R[2][0]
        +R[0][2]*R[1][0]*R[2][1]
        -R[0][2]*R[1][1]*R[2][0]
        -R[0][1]*R[1][0]*R[2][2]
        -R[0][0]*R[1][2]*R[2][1];
    if ( det < 0.0 )
    {
        R[0][0] = -R[0][0];
        R[1][0] = -R[1][0];
        R[2][0] = -R[2][0];
    }

    // extract angles from rotation axis = (cos(u)sin(v),sin(u)sin(v),cos(v))
    double axis[3];
    MatrixToAngleAxis(R,angle[2],axis);
    if ( -1+epsilon < axis[2] )
    {
        if ( axis[2] < 1-epsilon )
        {
            angle[0] = atan2(axis[1],axis[0]);
            angle[1] = acos(axis[2]);
        }
        else
        {
            angle[0] = 0;
            angle[1] = 0;
        }
    }
    else
    {
        angle[0] = 0;
        angle[1] = AGT_PI;
    }
}

Definition at line 134 of file AnalyticGeometryTool.cpp.

{
  if (instance_ == 0) {
    instance_ = new AnalyticGeometryTool;
  }
  return instance_;
}
int AnalyticGeometryTool::int_ln_ln ( double  p1[3],
double  v1[3],
double  p2[3],
double  v2[3],
double  int_pnt1[3],
double  int_pnt2[3] 
)

This function finds the intersection of two lines. If the lines cross, it finds the one unique point. If the lines do not cross (but are non- parallel), it finds the nearest intersection points (a line drawn from each of these intersections would be perpendicular to both lines). Returns number of intersections found: 0 intersections found, lines are parallel 1 intersection found, lines physically intersect 2 intersections found, lines do not intersect but nearest intersections found

Definition at line 1410 of file AnalyticGeometryTool.cpp.

{   
   double norm[3];
   double pln1_norm[3];
   double pln2_norm[3];
   
   // Cross the two vectors to get a normal vector 
   cross_vec(v1, v2, norm);
   
   if (dbl_eq(mag_vec(norm), 0.0))
      return 0;
   
   // Cross v2 & normal to get normal to plane 2 
   cross_vec(v2, norm, pln2_norm);
   
   // Find intersection of pln2_norm and vector 1 - this is int_pnt1 
   int_ln_pln(p1, v1, p2, pln2_norm, int_pnt1);
   
   // Cross v1 & normal to get normal to plane 1 
   cross_vec(v1, norm, pln1_norm);
   
   // Find intersection of pln2_norm and vector 1 - this is int_pnt2 
   int_ln_pln(p2, v2, p1, pln1_norm, int_pnt2);

   // Check to see if the intersection points are the same & return 
   if (pnt_eq(int_pnt1, int_pnt2)) 
      return 1;
   else
      return 2;   
}
CubitStatus AnalyticGeometryTool::int_ln_pln ( double  ln_orig[3],
double  ln_vec[3],
double  pln_orig[3],
double  pln_norm[3],
double  int_pnt[3] 
)

This routine calculates the intersection point of a line and a plane. Returns CUBIT_FAILURE if no intersection exists (line is parallel to the plane).

Definition at line 1381 of file AnalyticGeometryTool.cpp.

{   
   double denom;
   double t;
   
   // Set parametric eqns of line equal to parametric eqn of plane & solve
   // for t 
   denom = pln_norm[0]*ln_vec[0] + pln_norm[1]*ln_vec[1] + 
           pln_norm[2]*ln_vec[2];
   
   if (dbl_eq(denom, 0.0))
      return CUBIT_FAILURE;
   
   t = (pln_norm[0]*(pln_orig[0]-ln_orig[0]) + 
      pln_norm[1]*(pln_orig[1]-ln_orig[1]) + 
      pln_norm[2]*(pln_orig[2]-ln_orig[2]))/denom;
   
   // Substitute t back into equations of line to get xyz 
   int_pnt[0] = ln_orig[0] + ln_vec[0]*t;
   int_pnt[1] = ln_orig[1] + ln_vec[1]*t;
   int_pnt[2] = ln_orig[2] + ln_vec[2]*t;
      
   return CUBIT_SUCCESS;   
}
CubitStatus AnalyticGeometryTool::int_pln_pln ( double  pln_1_orig[3],
double  pln_1_norm[3],
double  pln_2_orig[3],
double  pln_2_norm[3],
double  ln_orig[3],
double  ln_vec[3] 
)

This routine finds intersection of two planes. This intersection results in a line. Returns CUBIT_FAILURE if planes are parallel.

int AnalyticGeometryTool::int_pnt_ln ( double  pnt[3],
double  ln_orig[3],
double  ln_vec[3],
double  int_pnt[3] 
)

This routine finds the perpendicular intersection of a point & line. This point will lie on the line. Returns 0 if point is not on the line, otherwise 1.

int AnalyticGeometryTool::int_pnt_pln ( double  pnt[3],
double  pln_orig[3],
double  pln_norm[3],
double  pln_int[3] 
)

This routine determines the perpendicular intersection of a point and a plane. This is the perpendicular intersection point on the plane. Returns 0 if point is not on the plane, otherwise 1.

Definition at line 1444 of file AnalyticGeometryTool.cpp.

{
   // Calculate line plane intersection w/plane normal as line vector 
   int_ln_pln(pnt, pln_norm, pln_orig, pln_norm, pln_int);
   
   // Check to see if point is on the plane 
   if (pnt_eq(pln_int, pnt))
      return 1;
   else
      return 0; 
}

Definition at line 6087 of file AnalyticGeometryTool.cpp.

{
    // build edges of T0
    AgtLine line[3];
    TriangleLines(T0,line);

    // initial list is copy of T1 (since triangle may be deleted)
    AgtTriList* list = new AgtTriList;
    list->tri = new Triangle;
    memcpy(list->tri,T1,sizeof(Triangle));
    list->next = 0;

    // process subtriangles of T1 against lines of T0
    for (int i = 0; i < 3; i++)
    {
        AgtTriList* save = 0;
        while ( list )
        {
            // get head of list
            Triangle* T = list->tri;
            
            AgtTriList* sad = SplitAndDecompose(T,&line[i]);
            if ( sad )
            {
                // search for end of list
                AgtTriList* end = sad;
                while ( end->next )
                    end = end->next;

                // attach decomposition to front of save-list
                end->next = save;
                save = sad;
            }

            // remove head of list
            AgtTriList* tmp = list;
            list = list->next;

            if( sad == NULL )
              delete tmp->tri;
            delete tmp;
        }
        list = save;
    }

    return list;
}
CubitStatus AnalyticGeometryTool::inv_mtx_adj ( double  mtx[3][3],
double  inv_mtx[3][3] 
)

This routine inverts a 3x3 matrix using the adjoint method. If NULL is sent in for first matrix, the second matrix is set to the identity matrix. If the same memory address is passed in for both incoming and outgoing matrices, the matrix is inverted in place. Returns CUBIT_FAILURE if no inverse exists.

Definition at line 808 of file AnalyticGeometryTool.cpp.

{
   int i,i1,i2,j,j1,j2;
   double work_mtx[3][3];
   double determinant;
   
   // Check for null input 
   if (!mtx) {
      copy_mtx(AGT_IDENTITY_MTX_3X3, inv_mtx);
      return CUBIT_SUCCESS;
   }      
   
   // Calculate determinant 
   determinant = det_mtx(mtx);
   
   // Check for singularity 
   if (dbl_eq(determinant,0.0))
      return CUBIT_FAILURE;
   
   // Get work matrix (allow inverting in place) 
   copy_mtx(mtx, work_mtx);
   
   // Inverse is adjoint matrix divided by determinant 
   for (i=1; i<4; i++) {
      
      i1 = (i % 3) + 1;  i2 = (i1 % 3) + 1;
      
      for (j=1; j<4; j++) {
         
         j1 = (j % 3) + 1;  j2 = (j1 % 3) + 1;
         
         inv_mtx[j-1][i-1] = (work_mtx[i1-1][j1-1]*work_mtx[i2-1][j2-1] - 
            work_mtx[i1-1][j2-1]*work_mtx[i2-1][j1-1]) 
            / determinant;
         
      }
   }
   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::inv_trans_mtx ( double  transf[4][4],
double  inv_transf[4][4] 
)

This routine inverts a 4x4 transformation matrix. If NULL is sent in for first matrix, the second matrix is set to the identity matrix. If the same memory address is passed in for both incoming and outgoing matrices, the matrix is inverted in place. Uses LU decomposition.

Definition at line 849 of file AnalyticGeometryTool.cpp.

{
   double scale_sq;
   double inv_sq_scale;
   double tmp_transf[4][4]; // For temporary storage of incoming matrix
   double *p_tmp_transf = NULL;

   // If input transform is 0 set output to identity matrix
   if (!transf) {
      copy_mtx( AGT_IDENTITY_MTX_4X4, inv_transf );
      return CUBIT_SUCCESS;
   }
   
   // Obtain the matrix scale
   scale_sq = transf[0][0]*transf[0][0] + transf[0][1]*transf[0][1] +
              transf[0][2]*transf[0][2];
   
   // Check for singular matrix
   if (scale_sq < (.000000001 * .000000001))
      return CUBIT_FAILURE;
   
   // Need the inverse scale squared
   inv_sq_scale = 1.0 / scale_sq;   
   
   // Check to see if inverting in place
   if (transf == inv_transf) {
      copy_mtx( transf, tmp_transf );
      p_tmp_transf = (double *)tmp_transf;
   }
   else
      p_tmp_transf = (double *)transf;
   
   // The X vector
   inv_transf[0][0] = p_tmp_transf[0] * inv_sq_scale;
   inv_transf[1][0] = p_tmp_transf[1] * inv_sq_scale;
   inv_transf[2][0] = p_tmp_transf[2] * inv_sq_scale;
   
   // The Y vector
   inv_transf[0][1] = p_tmp_transf[4] * inv_sq_scale;
   inv_transf[1][1] = p_tmp_transf[5] * inv_sq_scale;
   inv_transf[2][1] = p_tmp_transf[6] * inv_sq_scale;
   
   // The Z vector 
   inv_transf[0][2] = p_tmp_transf[8] * inv_sq_scale;
   inv_transf[1][2] = p_tmp_transf[9] * inv_sq_scale;
   inv_transf[2][2] = p_tmp_transf[10] * inv_sq_scale;
   
   // Column 4 
   inv_transf[0][3] = 0.0;
   inv_transf[1][3] = 0.0;
   inv_transf[2][3] = 0.0;
   
   // The X origin 
   inv_transf[3][0] = -inv_sq_scale * (  p_tmp_transf[0] * p_tmp_transf[12]
                       + p_tmp_transf[1] * p_tmp_transf[13]
                       + p_tmp_transf[2] * p_tmp_transf[14]);
   
   // The Y origin 
   inv_transf[3][1] = -inv_sq_scale * (  p_tmp_transf[4] * p_tmp_transf[12]
                       + p_tmp_transf[5] * p_tmp_transf[13]
                       + p_tmp_transf[6] * p_tmp_transf[14]);
   
   // The Z origin 
   inv_transf[3][2] = -inv_sq_scale * (  p_tmp_transf[8] * p_tmp_transf[12]
                       + p_tmp_transf[9] * p_tmp_transf[13]
                       + p_tmp_transf[10] * p_tmp_transf[14]);
   
   // This is always one 
   inv_transf[3][3] = 1.0;
   
   return CUBIT_SUCCESS;
}
CubitBoolean AnalyticGeometryTool::is_ln_on_pln ( double  ln_orig[3],
double  ln_vec[3],
double  pln_orig[3],
double  pln_norm[3] 
)

This routine determines if a specified line (defined by a point and a vector) is in a given plane (defined by the two vectors in the plane and the normal to that plane).

Definition at line 1700 of file AnalyticGeometryTool.cpp.

{     
   
   // Check to see if line origin is on the plane.  If so, check to see if
   // line vector is perpendicular to plane normal.
   
   if (is_pnt_on_pln(ln_orig, pln_orig, pln_norm) &&
       is_vec_perp(ln_vec, pln_norm))      
      return CUBIT_TRUE;      
   else      
      return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_pln_on_pln ( double  pln_orig1[3],
double  pln_norm1[3],
double  pln_orig2[3],
double  pln_norm2[3] 
)

This routine checks to see if two infinite planes are coincident.

Definition at line 1718 of file AnalyticGeometryTool.cpp.

{  
   // If 1st plane origin is on the 2nd plane and the normals are
   // parallel they are coincident
   if( is_vec_par( pln_norm1, pln_norm2 ) &&
      is_pnt_on_pln( pln_orig1, pln_orig2, pln_norm2 ) )
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_pnt_on_ln ( double  pnt[3],
double  ln_orig[3],
double  ln_vec[3] 
)

This routined determines if a specified point is on a given infinite line defined by a point and a vector.

Definition at line 1498 of file AnalyticGeometryTool.cpp.

{  
   double vec[3];
   
   // Compare pnt and line origin
   if (pnt_eq(pnt, ln_orig))
      return CUBIT_TRUE;   
   
   // Get a vector from line origin to the point      
   get_vec(ln_orig, pnt, vec);
   
   // If this vector is parallel with line vector, point is on the line
   if (is_vec_par(vec, ln_vec))
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_pnt_on_ln_seg ( double  pnt1[3],
double  end1[3],
double  end2[3] 
)

This routine determines if a specified point is on a given line segment defined by two endpoints. The point is on the line only if it lies on the line between or on the two endpoints.

Definition at line 1518 of file AnalyticGeometryTool.cpp.

{ 
//   METHOD:
//    o Use parametric equations of line
//    
//         x = x1 + t(x2 - x1)
//         y = y1 + t(y2 - y1)
//         z = z1 + t(z2 - z1)
//   
//    o Note:  two other method's were considered: 
//       1) Comparing sum of distance of point to both end points to the
//          line length.
//       2) Checking to see if area of a triangle with the vertices is zero
//       
//     Using parametric equations is more efficient in many cases.
  double t1 = 0.0,
    t2 = 0.0,
    t3 = 0.0,
    neg_range,
    pos_range;
  
  unsigned short flg1 = 0,
    flg2 = 0,
    flg3 = 0;
  
  neg_range = 0.0 - agtEpsilon;
  pos_range = 1.0 + agtEpsilon;
  
  if (fabs(end2[0] - end1[0]) < agtEpsilon)
  {
    if (fabs(pnt[0] - end1[0]) < agtEpsilon)
      flg1 = 1;
    else
      return CUBIT_FALSE;
  }
  else
  {
    t1 = (pnt[0] - end1[0])/(end2[0] - end1[0]);
    
    if (t1<neg_range || t1>pos_range)
      return CUBIT_FALSE;
  }
  
  if (fabs(end2[1] - end1[1]) < agtEpsilon)
  {
    if (fabs(pnt[1] - end1[1]) < agtEpsilon)
      flg2 = 1;
    else
      return CUBIT_FALSE;
  }
  else
  {
    t2 = (pnt[1] - end1[1])/(end2[1] - end1[1]);
    
    if (t2<neg_range || t2>pos_range)
      return CUBIT_FALSE;
  }
   
  if (fabs(end2[2] - end1[2]) < agtEpsilon)
  {
    if (fabs(pnt[2] - end1[2]) < agtEpsilon)
      flg3 = 1;
    else
      return CUBIT_FALSE;
  }
  else
  {
    t3 = (pnt[2] - end1[2])/(end2[2] - end1[2]);
    
    if (t3<neg_range || t3>pos_range)
      return CUBIT_FALSE;
  }
  
    // If any 2 flags are 1, point is on the line,
    // otherwise, check remaining T's for equality
  
  if (flg1)
  {
      // Here, flg1 = 1
    
    if (flg2)
    {
        // Here, flg1 = 1
        // Here, flg2 = 1
      return CUBIT_TRUE;
    }
    else
    {
        // Here, flg1 = 1
        // Here, flg2 = 0
      
      if (flg3)
          // Here, flg1 = 1
          // Here, flg2 = 0
          // Here, flg3 = 1
        return CUBIT_TRUE;
      else
      {
          // Here, flg1 = 1
          // Here, flg2 = 0
          // Here, flg3 = 0
        if (dbl_eq(t2, t3))
          return CUBIT_TRUE;
        else
          return CUBIT_FALSE;
      }
    }
  }
  else
  {
      // Here, flg1 = 0
    if (flg2)
    {
        // Here, flg1 = 0
        // Here, flg2 = 1
      if (flg3)
        return CUBIT_TRUE;
          // Here, flg1 = 0
          // Here, flg2 = 1
          // Here, flg3 = 1
      else
      {
          // Here, flg1 = 0
          // Here, flg2 = 1
          // Here, flg3 = 0
        if (dbl_eq(t1, t3))
          return CUBIT_TRUE;
        else
          return CUBIT_FALSE;
      }
    }
    else
    {
        // Here, flg1 = 0
        // Here, flg2 = 0
      if (flg3)
      {
          // Here, flg1 = 0
          // Here, flg2 = 0
          // Here, flg3 = 1
        if (dbl_eq(t1, t2))
          return CUBIT_TRUE;
        else
          return CUBIT_FALSE;
      }
      else
      {
          // Here, flg1 = 0
          // Here, flg2 = 0
          // Here, flg3 = 0
        if (dbl_eq(t1, t2) && dbl_eq(t1, t3))
          return CUBIT_TRUE;
        else
          return CUBIT_FALSE;
      }
    }
  }
  
    // This would be a programmer's error if we got to this point
//  return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_pnt_on_pln ( double  pnt[3],
double  pln_orig[3],
double  pln_norm[3] 
)

This routined determines if a specified point is on a given plane which is defined by an origin and three vectors.

Definition at line 1682 of file AnalyticGeometryTool.cpp.

{  
   double result;
   
   // See if point satisfies parametric equation of plane
   
   result = pln_norm[0] * (pnt[0] - pln_orig[0]) +
            pln_norm[1] * (pnt[1] - pln_orig[1]) +
            pln_norm[2] * (pnt[2] - pln_orig[2]);
   
   if (dbl_eq(result, 0.0))
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_vec_par ( double  vec_1[3],
double  vec_2[3] 
)

This routine checks to see if two vectors are parallel. Two vectors are considered parallel if they are pointing in the same direction or opposite directions.

Definition at line 1464 of file AnalyticGeometryTool.cpp.

{   
   double cross[3];
   
   // Get cross product & see if its magnitude is zero
   cross_vec(vec_1, vec_2, cross);
   
   if (dbl_eq(mag_vec(cross), 0.0))
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;
}
CubitBoolean AnalyticGeometryTool::is_vec_perp ( double  vec_1[3],
double  vec_2[3] 
)

This routine checks to see if two vectors are perpendicular. Two vectors are considered perpendicular if they are PI/2 radians to each other.

Definition at line 1478 of file AnalyticGeometryTool.cpp.

{  
   // Check angle between vectors
   if (dbl_eq(angle_vec_vec(vec_1, vec_2), AGT_PI_DIV_2))
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;   
}
CubitBoolean AnalyticGeometryTool::is_vecs_same_dir ( double  vec_1[3],
double  vec_2[3] 
)

This routine checks to see if two vectors point in the same direction. The vector magnitudes do not have to be equal for a successful return.

Definition at line 1487 of file AnalyticGeometryTool.cpp.

{         
   // Check to see if angle between vectors can be considered zero   
   if (dbl_eq(angle_vec_vec(vec_1, vec_2), 0.0))
      return CUBIT_TRUE;
   else
      return CUBIT_FALSE;   
}
double AnalyticGeometryTool::Length ( const Point2 p) [private]
double AnalyticGeometryTool::Length ( const Point3 p) [private]
double AnalyticGeometryTool::mag_vec ( double  vec[3])

Finds the magnitude of a vector: magnitude = sqrt (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);

Definition at line 1196 of file AnalyticGeometryTool.cpp.

{
   return (sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]));
}
void AnalyticGeometryTool::MatrixToAngleAxis ( double **  R,
double &  angle,
double  axis[3] 
) [private]

Definition at line 2722 of file AnalyticGeometryTool.cpp.

{
    // Let (x,y,z) be the unit-length axis and let A be an angle of rotation.
    // The rotation matrix is R = I + sin(A)*P + (1-cos(A))*P^2 where
    // I is the identity and
    //
    //       +-        -+
    //   P = |  0 +z -y |
    //       | -z  0 +x |
    //       | +y -x  0 |
    //       +-        -+
    //
    // Some algebra will show that
    //
    //   cos(A) = (trace(R)-1)/2  and  R - R^t = 2*sin(A)*P
    //
    // In the event that A = pi, R-R^t = 0 which prevents us from extracting
    // the axis through P.  Instead note that R = I+2*P^2 when A = pi, so
    // P^2 = (R-I)/2.  The diagonal entries of P^2 are x^2-1, y^2-1, and z^2-1.
    // We can solve these for axis (x,y,z).  Because the angle is pi, it does
    // not matter which sign you choose on the square roots.

    double trace = R[0][0]+R[1][1]+R[2][2];
    double cs = 0.5*(trace-1.0);
    if ( -1 < cs )
    {
        if ( cs < 1 )
            angle = acos(cs);
        else
            angle = 0;
    }
    else
    {
        angle = AGT_PI;
    }

    axis[0] = R[1][2]-R[2][1];
    axis[1] = R[2][0]-R[0][2];
    axis[2] = R[0][1]-R[1][0];
    double length = sqrt(axis[0]*axis[0]+axis[1]*axis[1]+axis[2]*axis[2]);
    const double epsilon = 1e-06;
    if ( length > epsilon )
    {
        axis[0] /= length;
        axis[1] /= length;
        axis[2] /= length;
    }
    else  // angle is 0 or pi
    {
        if ( angle > 1.0 )  // any number strictly between 0 and pi works
        {
            // angle must be pi
            axis[0] = sqrt(0.5*(1.0+R[0][0]));
            axis[1] = sqrt(0.5*(1.0+R[1][1]));
            axis[2] = sqrt(0.5*(1.0+R[2][2]));

            // determine signs of axis components
            double test[3];
            test[0] = R[0][0]*axis[0]+R[0][1]*axis[1]+R[0][2]*axis[2]-axis[0];
            test[1] = R[1][0]*axis[0]+R[1][1]*axis[1]+R[1][2]*axis[2]-axis[1];
            test[2] = R[2][0]*axis[0]+R[2][1]*axis[1]+R[2][2]*axis[2]-axis[2];
            length = test[0]*test[0]+test[1]*test[1]+test[2]*test[2];
            if ( length < epsilon )
                return;

            axis[1] = -axis[1];
            test[0] = R[0][0]*axis[0]+R[0][1]*axis[1]+R[0][2]*axis[2]-axis[0];
            test[1] = R[1][0]*axis[0]+R[1][1]*axis[1]+R[1][2]*axis[2]-axis[1];
            test[2] = R[2][0]*axis[0]+R[2][1]*axis[1]+R[2][2]*axis[2]-axis[2];
            length = test[0]*test[0]+test[1]*test[1]+test[2]*test[2];
            if ( length < epsilon )
                return;

            axis[2] = -axis[2];
            test[0] = R[0][0]*axis[0]+R[0][1]*axis[1]+R[0][2]*axis[2]-axis[0];
            test[1] = R[1][0]*axis[0]+R[1][1]*axis[1]+R[1][2]*axis[2]-axis[1];
            test[2] = R[2][0]*axis[0]+R[2][1]*axis[1]+R[2][2]*axis[2]-axis[2];
            length = test[0]*test[0]+test[1]*test[1]+test[2]*test[2];
            if ( length < epsilon )
                return;

            axis[1] = -axis[1];
            test[0] = R[0][0]*axis[0]+R[0][1]*axis[1]+R[0][2]*axis[2]-axis[0];
            test[1] = R[1][0]*axis[0]+R[1][1]*axis[1]+R[1][2]*axis[2]-axis[1];
            test[2] = R[2][0]*axis[0]+R[2][1]*axis[1]+R[2][2]*axis[2]-axis[2];
            length = test[0]*test[0]+test[1]*test[1]+test[2]*test[2];
            if ( length < epsilon )
                return;
        }
        else
        {
            // Angle is zero, matrix is the identity, no unique axis, so
            // return (1,0,0) for as good a guess as any.
            axis[0] = 1.0;
            axis[1] = 0.0;
            axis[2] = 0.0;
        }
    }
}
CubitStatus AnalyticGeometryTool::min_cyl_box_int ( double  radius,
CubitVector axis,
CubitVector center,
CubitBox box,
int  extension_type,
double  extension,
CubitVector start,
CubitVector end,
int  num_tess_pnts = 2048 
)

Finds the start and end of a cylinder that just intersects the input box (defined by bottom-left and top-right corners, axis aligned with the cubit coordinate system). The resultant cylinder's start and end points can be extended out by a percentage of cylinder's length or absolute (making it longer in both directions). extension_type - 0=none, 1=percentage, 2=absolute The num_tess_pnts is the number of line points along the circle the function projects to the box to calculate the minimum intersection (more points could give a more accurate result for non-axis aligned cylinders).

Definition at line 2358 of file AnalyticGeometryTool.cpp.

{
   CubitVector box_min = box.minimum();
   CubitVector box_max = box.maximum();
   
   double box_min_pnt[3], box_max_pnt[3], axis_vec[3], center_pnt[3];
   box_min.get_xyz( box_min_pnt ); box_max.get_xyz( box_max_pnt );
   axis.get_xyz( axis_vec ); center.get_xyz( center_pnt );
   
   double start_pnt[3], end_pnt[3];
   
   if( min_cyl_box_int( radius, axis_vec, center_pnt, 
                        box_min_pnt, box_max_pnt, 
                        extension_type, extension, 
                        start_pnt, end_pnt, num_tess_pnts )
                        == CUBIT_FAILURE )
     return CUBIT_FAILURE;
   
   start.set( start_pnt ); end.set( end_pnt );
   
   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::min_cyl_box_int ( double  radius,
double  axis[3],
double  center[3],
double  box_min[3],
double  box_max[3],
int  extension_type,
double  extension,
double  start[3],
double  end[3],
int  num_tess_pnts = 2048 
)

Finds the start and end of a cylinder that just intersects the input box (defined by bottom-left and top-right corners, axis aligned with the cubit coordinate system). The resultant cylinder's start and end points can be extended out by a percentage of cylinder's length or absolute (making it longer in both directions). extension_type - 0=none, 1=percentage, 2=absolute The num_tess_pnts is the number of line points along the circle the function projects to the box to calculate the minimum intersection (more points could give a more accurate result for non-axis aligned cylinders).

Definition at line 2391 of file AnalyticGeometryTool.cpp.

{
  double cyl_z[3];
  unit_vec( axis, cyl_z );

  //PRINT_INFO( "Axis = %f, %f, %f\n", cyl_z[0], cyl_z[1], cyl_z[2] );
  //PRINT_INFO( "Center = %f, %f, %f\n", center[0], center[1], center[2] );

  // Find transformation matrix to take a point into cylinder's
  // coordinate system
  double cubit2cyl_mtx[4][4], cyl2cubit_mtx[4][4];
  double cyl_x[3], cyl_y[3];
  orth_vecs( cyl_z, cyl_x, cyl_y );
  vecs_to_mtx( cyl_x, cyl_y, cyl_z, center, cyl2cubit_mtx );
  inv_trans_mtx( cyl2cubit_mtx, cubit2cyl_mtx );

  // Setup the circle
  double vec_1[3], vec_2[3];
  orth_vecs( cyl_z, vec_1, vec_2 );
  AGT_3D_Arc arc;
  setup_arc( vec_1, vec_2, center, 0.0, 2.0 * CUBIT_PI, radius, arc );

  // Setup the planes of the box
  double pln_norm[6][3], pln_orig[6][3];
  // Front
  pln_orig[0][0] = 0.0; pln_orig[0][1] = 0.0; pln_orig[0][2] = box_max[2];
  pln_norm[0][0] = 0.0; pln_norm[0][1] = 0.0; pln_norm[0][2] = 1.0;
  // Left
  pln_orig[1][0] = box_min[0]; pln_orig[1][1] = 0.0; pln_orig[1][2] = 0.0;
  pln_norm[1][0] = -1.0; pln_norm[1][1] = 0.0; pln_norm[1][2] = 0.0;
  // Top
  pln_orig[2][0] = 0.0; pln_orig[2][1] = box_max[1]; pln_orig[2][2] = 0.0;
  pln_norm[2][0] = 0.0; pln_norm[2][1] = 1.0; pln_norm[2][2] = 0.0;
  // Right
  pln_orig[3][0] = box_max[0]; pln_orig[3][1] = 0.0; pln_orig[3][2] = 0.0;
  pln_norm[3][0] = 1.0; pln_norm[3][1] = 0.0; pln_norm[3][2] = 0.0;
  // Bottom
  pln_orig[4][0] = 0.0; pln_orig[4][1] = box_min[1]; pln_orig[4][2] = 0.0;
  pln_norm[4][0] = 0.0; pln_norm[4][1] = -1.0; pln_norm[4][2] = 0.0;
  // Back
  pln_orig[5][0] = 0.0; pln_orig[5][1] = 0.0; pln_orig[5][2] = box_min[2];
  pln_norm[5][0] = 0.0; pln_norm[5][1] = 0.0; pln_norm[5][2] = -1.0;

  double z; // Intersection along cylinder's axis
  double min_z = CUBIT_DBL_MAX, max_z = -CUBIT_DBL_MAX;
  
  double t = 0.0, dt;
  dt = 1.0/(double)num_tess_pnts;
  double pnt[3];
  double int_pnt[3];
  double box_tol = 1e-14;
  double box_min_0 = box_min[0]-box_tol;
  double box_min_1 = box_min[1]-box_tol;
  double box_min_2 = box_min[2]-box_tol;
  double box_max_0 = box_max[0]+box_tol;
  double box_max_1 = box_max[1]+box_tol;
  double box_max_2 = box_max[2]+box_tol;

  int i,j;
  for( i=0; i<num_tess_pnts; i++ )
  {
    get_arc_xyz( arc, t, pnt );

    for( j=0; j<6; j++ )
    {
      // Evaluate the intersection at this point
      if( int_ln_pln( pnt, cyl_z, pln_orig[j], pln_norm[j], int_pnt )
        == CUBIT_FAILURE )
        continue;

      // Throw-out if intersection not on physical box
      if( int_pnt[0] < box_min_0 || int_pnt[1] < box_min_1 ||
          int_pnt[2] < box_min_2 || int_pnt[0] > box_max_0 ||
          int_pnt[1] > box_max_1 || int_pnt[2] > box_max_2 ) 
        continue;

      // Find min/max cylinder z on box so far
      // z-distance (in cylinder coordinate system)
      z = cubit2cyl_mtx[0][2]*int_pnt[0] + cubit2cyl_mtx[1][2]*int_pnt[1] + 
        cubit2cyl_mtx[2][2]*int_pnt[2] + cubit2cyl_mtx[3][2];

      if( z < min_z ) min_z = z;
      if( z > max_z ) max_z = z;
      
    }

    t += dt;

  }

  // Check the 8 corners of the box - they are likely min/max's.
  double box_corners[8][3];
  get_box_corners( box_min, box_max, box_corners );
  for( i=0; i<8; i++ )
  {
    // Get the corner in the cylinder csys
    transform_pnt( cubit2cyl_mtx, box_corners[i], pnt );
    // If the pnt is within the circle's radius, check it's z-coord
    // (distance from center)
    if(  sqrt( pnt[0]*pnt[0] + pnt[1]*pnt[1] ) <= radius+box_tol )
    {
      if( pnt[2] < min_z ) min_z = pnt[2];
      if( pnt[2] > max_z ) max_z = pnt[2];
    }
  }

  if( min_z == CUBIT_DBL_MAX || max_z == -CUBIT_DBL_MAX )
  {
    PRINT_ERROR( "Unable to find cylinder/box intersection\n" );
    return CUBIT_FAILURE;
  }

  if( min_z == max_z )
  {
    PRINT_ERROR( "Unable to find cylinder/box intersection\n" );
    return CUBIT_FAILURE;
  }

  //PRINT_INFO( "Min dist = %f\n", min_z );
  //PRINT_INFO( "Max dist = %f\n", max_z );

  // Find the start and end of the cylinder
  next_pnt( center, cyl_z, min_z, start );
  next_pnt( center, cyl_z, max_z, end );

  //PRINT_INFO( "Start = %f, %f, %f\n", start[0], start[1], start[2] );
  //PRINT_INFO( "End = %f, %f, %f\n", end[0], end[1], end[2] );

  // Extend start and end, if necessary
  if( extension_type > 0 )
  {
    double ext_distance = 0.0;
    if( extension_type == 1 ) // percentage
    {
      double cyl_length = dist_pnt_pnt( start, end );
      ext_distance = extension/100.0 * cyl_length;
    }
    else
      ext_distance = extension;

    next_pnt( end, cyl_z, ext_distance, end );
    reverse_vec( cyl_z, cyl_z );
    next_pnt( start, cyl_z, ext_distance, start );
  }

  return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::min_pln_box_int_corners ( const CubitPlane plane,
const CubitBox box,
int  extension_type,
double  extension,
CubitVector p1,
CubitVector p2,
CubitVector p3,
CubitVector p4,
CubitBoolean  silent = CUBIT_FALSE 
)

Finds the 4 corner points of the input infinite plane that just intersects the input box (defined by bottom-left and top-right corners, axis aligned with the cubit coordinate system) - plane should have minimal area. The resultant plane's corner points can be extended out by a percentage of diagonal or absolute (making it bigger than the minimal area plane). extension_type - 0=none, 1=percentage, 2=absolute If silent is CUBIT_TRUE, no error messages are given. Note this returns points even if the plane doesn't physically intersect the given box... it does this by moving the plane to the center of the box, doing the intersection, then projecting the points back to the original plane.

Definition at line 1912 of file AnalyticGeometryTool.cpp.

{
   CubitVector box_min = box.minimum();
   CubitVector box_max = box.maximum();
   
   CubitVector plane_norm = plane.normal();
   
   double box_min_pnt[3], box_max_pnt[3], pln_norm[3];
   box_min.get_xyz( box_min_pnt ); box_max.get_xyz( box_max_pnt );
   plane_norm.get_xyz( pln_norm );
   
   double pnt1[3], pnt2[3], pnt3[3], pnt4[3];
   
   if( min_pln_box_int_corners( pln_norm, plane.coefficient(),
      box_min_pnt, box_max_pnt, 
      extension_type, extension, 
      pnt1, pnt2, pnt3, pnt4, silent ) == CUBIT_FAILURE )
         return CUBIT_FAILURE;
   
   p1.set( pnt1 ); p2.set( pnt2 ); p3.set( pnt3 ); p4.set( pnt4 );
   
   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::min_pln_box_int_corners ( CubitVector vec1,
CubitVector vec2,
CubitVector vec3,
CubitVector box_min,
CubitVector box_max,
int  extension_type,
double  extension,
CubitVector p1,
CubitVector p2,
CubitVector p3,
CubitVector p4,
CubitBoolean  silent = CUBIT_FALSE 
)

Finds the 4 corner points of the input infinite plane that just intersects the input box (defined by bottom-left and top-right corners, axis aligned with the cubit coordinate system) - plane should have minimal area. The resultant plane's corner points can be extended out by a percentage of diagonal or absolute (making it bigger than the minimal area plane). extension_type - 0=none, 1=percentage, 2=absolute If silent is CUBIT_TRUE, no error messages are given. Note this returns points even if the plane doesn't physically intersect the given box... it does this by moving the plane to the center of the box, doing the intersection, then projecting the points back to the original plane.

Definition at line 1943 of file AnalyticGeometryTool.cpp.

{
   CubitPlane plane;
   if( plane.mk_plane_with_points( vec1, vec2, vec3 ) == CUBIT_FAILURE )
      return CUBIT_FAILURE;

   CubitVector plane_norm = plane.normal();
   double coefficient = plane.coefficient();

   double plane_norm3[3];
   double box_min3[3];
   double box_max3[3];

   box_min.get_xyz( box_min3 );
   box_max.get_xyz( box_max3 );
   plane_norm.get_xyz( plane_norm3 );

   double p1_3[3], p2_3[3], p3_3[3], p4_3[3];
   p1.get_xyz( p1_3 );
   p2.get_xyz( p2_3 );
   p3.get_xyz( p3_3 );
   p4.get_xyz( p4_3 );
   
   if( min_pln_box_int_corners( plane_norm3, coefficient, box_min3, box_max3, 
      extension_type, extension, p1_3, p2_3, p3_3, p4_3, silent ) == CUBIT_FAILURE )
      return CUBIT_FAILURE;

   p1.set( p1_3 );
   p2.set( p2_3 );
   p3.set( p3_3 );
   p4.set( p4_3 );

   return CUBIT_SUCCESS;
}
CubitStatus AnalyticGeometryTool::min_pln_box_int_corners ( double  plane_norm[3],
double  plane_coeff,
double  box_min[3],
double  box_max[3],
int  extension_type,
double  extension,
double  p1[3],
double  p2[3],
double  p3[3],
double  p4[3],
CubitBoolean  silent = CUBIT_FALSE 
)

Finds the 4 corner points of the input infinite plane that just intersects the input box (defined by bottom-left and top-right corners, axis aligned with the cubit coordinate system) - plane should have minimal area. The resultant plane's corner points can be extended out by a percentage of diagonal or absolute (making it bigger than the minimal area plane). extension_type - 0=none, 1=percentage, 2=absolute If silent is CUBIT_TRUE, no error messages are given. Note this returns points even if the plane doesn't physically intersect the given box... it does this by moving the plane to the center of the box, doing the intersection, then projecting the points back to the original plane.

Definition at line 1988 of file AnalyticGeometryTool.cpp.

{
   int i;
   double cubit2pln_mtx[4][4],
      pln2cubit_mtx[4][4];
   double pln_orig[3];
   double box_tol = 1e-6;

   // Adjust bounding box if it is zero in any direction
   double abox_min[3], abox_max[3];
   copy_pnt( box_min, abox_min );
   copy_pnt( box_max, abox_max );

   double x_range = fabs(box_max[0]-box_min[0]);
   double y_range = fabs(box_max[1]-box_min[1]);
   double z_range = fabs(box_max[2]-box_min[2]);

   if( x_range < box_tol || y_range < box_tol || z_range < box_tol )
   {
     // Adjust zero length side(s) to maximum range
     double adj = CUBIT_MAX((CUBIT_MAX((x_range),(y_range))),(z_range))/2.0;

     // Arbitrarily expand box if it is zero in size
     if( adj < box_tol )
       adj = 5.0;  

     if( x_range < box_tol )
     {
       abox_min[0] -= adj;
       abox_max[0] += adj;
     }
     if( y_range < box_tol )
     {
       abox_min[1] -= adj;
       abox_max[1] += adj;
     }
     if( z_range < box_tol )
     {
       abox_min[2] -= adj;
       abox_max[2] += adj;
     }
   }

   // Get plane origin
   double A = pln_norm[0];
   double B = pln_norm[1];
   double C = pln_norm[2];
   double D = pln_coeff;
   get_pln_orig_norm( A, B, C, D, pln_orig );

   // Find intersections of edges with plane.  Add to unique array.  At most
   // there are 6 intersections...
   double int_array[6][3];
   int num_int = 0;
   num_int = get_plane_bbox_intersections( abox_min, abox_max, pln_orig, pln_norm,
     int_array );

   // Adjust so plane intercepts bounding box
   if( num_int < 3 )
   {
     // Move the plane to the center of the bounding box
     double box_cent[3];
     box_cent[0] = (abox_min[0]+abox_max[0])/2.0;
     box_cent[1] = (abox_min[1]+abox_max[1])/2.0;
     box_cent[2] = (abox_min[2]+abox_max[2])/2.0;

     // Get the intersections
     num_int = get_plane_bbox_intersections( abox_min, abox_max, box_cent, pln_norm, int_array );

     // Project the points back to the original plane
     switch (num_int)
     {
     case 6:
       int_pnt_pln( int_array[5], pln_orig, pln_norm, int_array[5] );
     case 5:
       int_pnt_pln( int_array[4], pln_orig, pln_norm, int_array[4] );
     case 4:
       int_pnt_pln( int_array[3], pln_orig, pln_norm, int_array[3] );
     case 3:
       int_pnt_pln( int_array[2], pln_orig, pln_norm, int_array[2] );
     case 2:
       int_pnt_pln( int_array[1], pln_orig, pln_norm, int_array[1] );
     case 1:
       int_pnt_pln( int_array[0], pln_orig, pln_norm, int_array[0] );
     }
   }

   if( num_int == 0 )
   {
     if( silent == CUBIT_FALSE )
       PRINT_ERROR( "Plane does not intersect the bounding box\n"
       "       Can't find 4 corners of plane\n" );
     return CUBIT_FAILURE;
   }
   if( num_int < 3 )
   {
     if( silent == CUBIT_FALSE )
       PRINT_ERROR( "Plane intersects the bounding box at only %d locations\n"
       "      Can't calculate 4 corners of plane\n", num_int );
      return CUBIT_FAILURE;
   }

   // Transform pnts to plane coordinate system
   double pln_x[3], pln_y[3];
   orth_vecs( pln_norm, pln_x, pln_y );
   vecs_to_mtx( pln_x, pln_y, pln_norm, pln_orig, pln2cubit_mtx );
   inv_trans_mtx( pln2cubit_mtx, cubit2pln_mtx );

   double int_arr_pln[6][3];
   for( i=0; i<num_int; i++ )
      transform_pnt( cubit2pln_mtx, int_array[i], int_arr_pln[i] );

   for(i=0; i<num_int; i++)
   {
     if( !dbl_eq( int_arr_pln[i][2], 0.0 ) )
     {
       if( silent == CUBIT_FALSE )
         PRINT_ERROR( "in AnalyticGeometryTool::min_box_pln_int_corners\n"
         "       Transform to plane wrong\n" );
        return CUBIT_FAILURE;
     }
   }

   // Place into format for mimimal box calculation
   Point2 pt[6];
   for ( i=0; i<num_int; i++ )
   {
      pt[i].x = int_arr_pln[i][0];
      pt[i].y = int_arr_pln[i][1];
   }
   
   // Find rectangle with minimal area to surround the points
   // (this is definitely overkill esp. for the simple cases.....)
   OBBox2 minimal = MinimalBox2( num_int, pt );
   
   // Strip out results
   double old_epsilon = set_epsilon( 1e-10 );
   double centroid[3];
   centroid[0] = minimal.center.x;
   centroid[1] = minimal.center.y;
   centroid[2] = 0.0;
   round_near_val( centroid[0] ); // Makes near -1, 0, 1 values -1, 0, 1
   round_near_val( centroid[1] );
   transform_pnt( pln2cubit_mtx, centroid, centroid );
   
   double x_axis[3];
   x_axis[0] = minimal.axis[0].x;
   x_axis[1] = minimal.axis[0].y;
   x_axis[2] = 0.0;
   round_near_val( x_axis[0] );
   round_near_val( x_axis[1] );
   transform_vec( pln2cubit_mtx, x_axis, x_axis );
   
   double y_axis[3];
   y_axis[0] = minimal.axis[1].x;
   y_axis[1] = minimal.axis[1].y;
   y_axis[2] = 0.0;
   round_near_val( y_axis[0] );
   round_near_val( y_axis[1] );
   transform_vec( pln2cubit_mtx, y_axis, y_axis );

   set_epsilon( old_epsilon );
   
   double dist_x;
   double dist_y;
   double extension_distance = 0.0;
   if( extension_type == 1 ) // Percentage (of 1/2 diagonal)
   {
      double diag_len = sqrt(  minimal.extent[0]*minimal.extent[0]
                             + minimal.extent[1]*minimal.extent[1] );
      extension_distance = diag_len*extension/100.0;
   }
   else if( extension_type == 2 ) // Absolute distance in x and y
      extension_distance = extension;

   dist_x = minimal.extent[0] + extension_distance;
   dist_y = minimal.extent[1] + extension_distance;
   
   next_pnt( centroid, x_axis, -dist_x, p1 );
   next_pnt( p1, y_axis, -dist_y, p1 );
   
   next_pnt( centroid, x_axis, -dist_x, p2 );
   next_pnt( p2, y_axis, dist_y, p2 );
   
   next_pnt( centroid, x_axis, dist_x, p3 );
   next_pnt( p3, y_axis, dist_y, p3 );
   
   next_pnt( centroid, x_axis, dist_x, p4 );
   next_pnt( p4, y_axis, -dist_y, p4 );
   
   return CUBIT_SUCCESS;
}
OBBox2 AnalyticGeometryTool::MinimalBox2 ( int  N,
Point2 pt 
) [private]

Definition at line 2607 of file AnalyticGeometryTool.cpp.

{
    OBBox2 box;

    // bracket a minimum for angles in [-pi,pi]
    double angle, area;
    int imin = 0;
    double areaMin = Area(N,pt,angleMin);

    int i;
    for (i = 1; i <= maxPartition; i++)
    {
        angle = angleMin+i*angleRange/maxPartition;
        area = Area(N,pt,angle);
        if ( area < areaMin )
        {
            imin = i;
            areaMin = area;
        }
    }

    double angle0 = angleMin+(imin-1)*angleRange/maxPartition;
    double area0 = Area(N,pt,angle0);
    double angle1 = angleMin+(imin+1)*angleRange/maxPartition;
    double area1 = Area(N,pt,angle1);
    
    // use inverse parabolic interpolation to find the minimum
    for (i = 0; i <= invInterp; i++)
    {
        double angleMid, areaMid;

        // test for convergence (do not change these parameters)
        const double epsilon = 1e-08, tol = 1e-04;
//        const double omtol = 1.0-tol;
        if ( fabs(angle1-angle0) <= 2*tol*fabs(angle)+epsilon )
            break;

        // compute vertex of interpolating parabola
        double dangle0 = angle0-angle, dangle1 = angle1-angle;
        double darea0 = area0-area, darea1 = area1-area;
        double temp0 = dangle0*darea1, temp1 = dangle1*darea0;
        double delta = temp1-temp0;
        if ( fabs(delta) < epsilon )
           break;

        angleMid = angle+0.5*(dangle1*temp1-dangle0*temp0)/(temp1-temp0);

        // update bracket
        if ( angleMid < angle )
        {
            areaMid = Area(N,pt,angleMid);
            if ( areaMid <= area )
            {
                angle1 = angle;
                area1 = area;
                angle = angleMid;
                area = areaMid;
            }
            else
            {
                angle0 = angleMid;
                area0 = areaMid;
            }
        }
        else if ( angleMid > angle )
        {
            areaMid = Area(N,pt,angleMid);
            if ( areaMid <= area )
            {
                angle0 = angle;
                area0 = area;
                angle = angleMid;
                area = areaMid;
            }
            else
            {
                angle1 = angleMid;
                area1 = areaMid;
            }
        }
        else
        {
            // bracket middle already vertex of parabola
            break;
        }
    }

    MinimalBoxForAngle(N,pt,angle,box);
    return box;
}
OBBox3 AnalyticGeometryTool::MinimalBox3 ( int  N,
Point3 pt 
) [private]

Definition at line 3263 of file AnalyticGeometryTool.cpp.

{
    // compute a good initial guess for an oriented bounding box
    double angle[3];
    InitialGuess(N,pt,angle);
    double oldVolume = Volume(N,pt,angle);
 
  //if initial guess gives angles that are zero, the rest of this code 
  //can give bounding box that isn't tight.  This probably 
  //isn't the best fix but it works.  CDE  4-20-2009
  if( angle[0] == 0 )
    angle[0] = 0.5;
  if( angle[1] == 0 )
    angle[1] = 0.5;
  if( angle[2] == 0 )
    angle[2] = 0.5;

    double saveAngle[3] = { angle[0], angle[1], angle[2] };

    // Powell's direction set method
    double U[3][3], volume;
    const int maxiters = 3*32;
    for (int iter = 0; iter < maxiters; iter++)
    {
        // reset directions to avoid linear dependence degeneration
        if ( iter % 3 == 0 )
        {
            U[0][0] = 1.0;  U[0][1] = 0.0;  U[0][2] = 0.0;
            U[1][0] = 0.0;  U[1][1] = 1.0;  U[1][2] = 0.0;
            U[2][0] = 0.0;  U[2][1] = 0.0;  U[2][2] = 1.0;
        }

        // find minima in specified directions
        for (int d = 0; d < 3; d++)
            volume = MinimizeOnInterval(N,pt,angle,U[d]);

        // estimate a conjugate direction
        double conj[3] =
        {
            angle[0]-saveAngle[0],
            angle[1]-saveAngle[1],
            angle[2]-saveAngle[2]
        };
        double length = sqrt(conj[0]*conj[0]+conj[1]*conj[1]+conj[2]*conj[2]);
        if ( length >= 1e-06 )
        {
            double invLen = 1.0/length;
            conj[0] *= invLen;
            conj[1] *= invLen;
            conj[2] *= invLen;
            
            // minimize in conjugate direction
            volume = MinimizeOnInterval(N,pt,angle,conj);
        }
        else
        {
            // Possible local, but not global, minimum.  Search nearby for
            // a smaller volume.
            volume = MinimizeOnLattice(N,pt,angle,2,0.0001);
            volume = MinimizeOnLattice(N,pt,angle,2,0.0010);
            volume = MinimizeOnLattice(N,pt,angle,2,0.0100);
            volume = MinimizeOnLattice(N,pt,angle,2,0.1000);
        }

        // test for convergence
        const double epsilon = 1e-04;
        double diff = fabs(volume-oldVolume);
        if ( diff <= epsilon )
        {
            // Possible local, but not global, minimum.  Search nearby for
            // a smaller volume.
            volume = MinimizeOnLattice(N,pt,angle,2,0.0001);
            volume = MinimizeOnLattice(N,pt,angle,2,0.0010);
            volume = MinimizeOnLattice(N,pt,angle,2,0.0100);
            volume = MinimizeOnLattice(N,pt,angle,2,0.1000);
            diff = fabs(volume-oldVolume);
            if ( diff <= epsilon )
                break;
        }

        // cycle the directions and add conjugate direction to list
        U[0][0] = U[1][0];  U[0][1] = U[1][1];  U[0][2] = U[1][2];
        U[1][0] = U[2][0];  U[1][1] = U[2][1];  U[1][2] = U[2][2];
        U[2][0] = conj[0];  U[2][1] = conj[1];  U[2][2] = conj[2];

        // set parameters for next pass
        oldVolume = volume;
        saveAngle[0] = angle[0];
        saveAngle[1] = angle[1];
        saveAngle[2] = angle[2];
    }

    OBBox3 box;
    MinimalBoxForAngles(N,pt,angle,box);
    return box;
}
void AnalyticGeometryTool::MinimalBoxForAngle ( int  N,
Point2 pt,
double  angle,
OBBox2 box 
) [private]

Definition at line 2572 of file AnalyticGeometryTool.cpp.

{
    double cs = cos(angle), sn = sin(angle);


    double umin = +cs*pt[0].x+sn*pt[0].y, umax = umin;
    double vmin = -sn*pt[0].x+cs*pt[0].y, vmax = vmin;
    for (int i = 1; i < N; i++)
    {
        double u = +cs*pt[i].x+sn*pt[i].y;
        if ( u < umin )
            umin = u;
        else if ( u > umax )
            umax = u;

        double v = -sn*pt[i].x+cs*pt[i].y;
        if ( v < vmin )
            vmin = v;
        else if ( v > vmax )
            vmax = v;
    }

    double umid = 0.5*(umax+umin);
    double vmid = 0.5*(vmax+vmin);
    box.center.x = umid*cs-vmid*sn;
    box.center.y = umid*sn+vmid*cs;
    box.axis[0].x = cs;
    box.axis[0].y = sn;
    box.axis[1].x = -sn;
    box.axis[1].y = cs;
    box.extent[0] = 0.5*(umax-umin);
    box.extent[1] = 0.5*(vmax-vmin);
}
void AnalyticGeometryTool::MinimalBoxForAngles ( int  N,
Point3 pt,
double  angle[3],
OBBox3 box 
) [private]

Definition at line 2891 of file AnalyticGeometryTool.cpp.

{
    double cs0 = cos(angle[0]), sn0 = sin(angle[0]);
    double cs1 = cos(angle[1]), sn1 = sin(angle[1]);
    double axis[3] = { cs0*sn1, sn0*sn1, cs1 };
    double rot[3][3];
    AngleAxisToMatrix(angle[2],axis,rot);

    double min[3] =
    {
        rot[0][0]*pt[0].x+rot[1][0]*pt[0].y+rot[2][0]*pt[0].z,
        rot[0][1]*pt[0].x+rot[1][1]*pt[0].y+rot[2][1]*pt[0].z,
        rot[0][2]*pt[0].x+rot[1][2]*pt[0].y+rot[2][2]*pt[0].z
    };

    double max[3] = { min[0], min[1], min[2] };

    for (int i = 1; i < N; i++)
    {
        double test[3] =
        {
            rot[0][0]*pt[i].x+rot[1][0]*pt[i].y+rot[2][0]*pt[i].z,
            rot[0][1]*pt[i].x+rot[1][1]*pt[i].y+rot[2][1]*pt[i].z,
            rot[0][2]*pt[i].x+rot[1][2]*pt[i].y+rot[2][2]*pt[i].z
        };

        if ( test[0] < min[0] )
            min[0] = test[0];
        else if ( test[0] > max[0] )
            max[0] = test[0];

        if ( test[1] < min[1] )
            min[1] = test[1];
        else if ( test[1] > max[1] )
            max[1] = test[1];

        if ( test[2] < min[2] )
            min[2] = test[2];
        else if ( test[2] > max[2] )
            max[2] = test[2];
    }

    double mid[3] =
    {
        0.5*(max[0]+min[0]), 0.5*(max[1]+min[1]), 0.5*(max[2]+min[2])
    };

    box.center.x = mid[0]*rot[0][0]+mid[1]*rot[0][1]+mid[2]*rot[0][2];
    box.center.y = mid[0]*rot[1][0]+mid[1]*rot[1][1]+mid[2]*rot[1][2];
    box.center.z = mid[0]*rot[2][0]+mid[1]*rot[2][1]+mid[2]*rot[2][2];
    box.axis[0].x = rot[0][0];
    box.axis[0].y = rot[1][0];
    box.axis[0].z = rot[2][0];
    box.axis[1].x = rot[0][1];
    box.axis[1].y = rot[1][1];
    box.axis[1].z = rot[2][1];
    box.axis[2].x = rot[0][2];
    box.axis[2].y = rot[1][2];
    box.axis[2].z = rot[2][2];
    box.extent[0] = 0.5*(max[0]-min[0]);
    box.extent[1] = 0.5*(max[1]-min[1]);
    box.extent[2] = 0.5*(max[2]-min[2]);
}
double AnalyticGeometryTool::MinimizeOnInterval ( int  N,
Point3 pt,
double  A[3],
double  D[3] 
) [private]

Definition at line 3015 of file AnalyticGeometryTool.cpp.

{
    // compute intersection of line A+t*D with domain of function
    double tmin, tmax;
    GetInterval(A,D,tmin,tmax);
    double tran = tmax-tmin;
    double angle[3];

    if( tran == 0.0 ) // Added by SRS 10-6-2000
       tran = 1.0;

    // bracket a minimum for angles in [A+tmin*D,A+tmax*D]
    double t = 0.0;
    double volumeMin = Volume(N,pt,A);
    double volume;

    const int max_partition = 64;
    int i, imin;
    for (i = 0, imin = -1; i <= max_partition; i++)
    {
        t = tmin+i*tran/max_partition;
        Combine(angle,A,t,D);

        volume = Volume(N,pt,angle);
        if ( volume < volumeMin )
        {
            imin = i;
            volumeMin = volume;
        }
    }

    if ( imin != -1 )
    {
        t = tmin+imin*tran/max_partition;
    }
    else
    {
        t = 0.0;

        // interval in which t=0 lies
        imin = int(-tmin*max_partition/tran+0.5);
    }
    volume = volumeMin;

    double t0 = tmin+(imin-1)*tran/max_partition;
    Combine(angle,A,t0,D);
    double volume0 = Volume(N,pt,angle);

    double t1 = tmin+(imin+1)*tran/max_partition;
    Combine(angle,A,t1,D);
    double volume1 = Volume(N,pt,angle);
    
    // use inverse parabolic interpolation to find the minimum
    const int inv_interp = 64;
    for (i = 0; i <= inv_interp; i++)
    {
        double tMid, volumeMid;

        // test for convergence (do not change these parameters)
        const double epsilon = 1e-08, tol = 1e-04;
//        const double omtol = 1.0-tol;
        if ( fabs(t1-t0) <= 2*tol*fabs(t)+epsilon )
            break;

        // compute vertex of interpolating parabola
        double dt0 = t0-t, dt1 = t1-t;
        double dvolume0 = volume0-volume, dvolume1 = volume1-volume;
        double temp0 = dt0*dvolume1, temp1 = dt1*dvolume0;
        double delta = temp1-temp0;
        if ( fabs(delta) < epsilon )
           break;

        tMid = t+0.5*(dt1*temp1-dt0*temp0)/(temp1-temp0);

        // update bracket
        if ( tMid < t )
        {
            Combine(angle,A,tMid,D);
            volumeMid = Volume(N,pt,angle);
            if ( volumeMid <= volume )
            {
                t1 = t;
                volume1 = volume;
                t = tMid;
                volume = volumeMid;
            }
            else
            {
                t0 = tMid;
                volume0 = volumeMid;
            }
        }
        else if ( tMid > t )
        {
            Combine(angle,A,tMid,D);
            volumeMid = Volume(N,pt,angle);
            if ( volumeMid <= volume )
            {
                t0 = t;
                volume0 = volume;
                t = tMid;
                volume = volumeMid;
            }
            else
            {
                t1 = tMid;
                volume1 = volumeMid;
            }
        }
        else
        {
            // bracket middle already vertex of parabola
            break;
        }
    }

    Combine(A,A,t,D);
    return volume;
}
double AnalyticGeometryTool::MinimizeOnLattice ( int  N,
Point3 pt,
double  A[3],
int  layers,
double  thickness 
) [private]

Definition at line 3135 of file AnalyticGeometryTool.cpp.

{
    int xmin = 0, ymin = 0, zmin = 0;
    double volume = Volume(N,pt,A);

    double angle[3];
    for (int z = -layers; z <= layers; z++)
    {
        angle[2] = A[2]+thickness*z/layers;
        for (int y = -layers; y <= layers; y++)
        {
            angle[1] = A[1]+thickness*y/layers;
            for (int x = -layers; x <= layers; x++)
            {
                angle[0] = A[0]+thickness*x/layers;

                double v = Volume(N,pt,angle);
                if ( v < volume )
                {
                    xmin = x;
                    ymin = y;
                    zmin = z;
                    volume = v;
                }
            }
        }
    }

    A[0] += thickness*xmin/layers;
    A[1] += thickness*ymin/layers;
    A[2] += thickness*zmin/layers;

    return volume;
}
double AnalyticGeometryTool::MinLineSegmentLineSegment ( const Line3 seg0,
const Line3 seg1,
double &  s,
double &  t 
) [private]

Definition at line 4584 of file AnalyticGeometryTool.cpp.

{
    Point3 diff = seg0.b - seg1.b;
    double A = Dot(seg0.m,seg0.m);
    double B = -Dot(seg0.m,seg1.m);
    double C = Dot(seg1.m,seg1.m);
    double D = Dot(seg0.m,diff);
    double E;  // -Dot(seg1.m,diff), defer until needed
    double F = Dot(diff,diff);
    double det = Abs(A*C-B*B);  // A*C-B*B = |Cross(M0,M1)|^2 >= 0

    double tmp;

    if ( det >= par_tolerance )
    {
        // line segments are not parallel
        E = -Dot(seg1.m,diff);
        s = B*E-C*D;
        t = B*D-A*E;
        
        if ( s >= 0 )
        {
            if ( s <= det )
            {
                if ( t >= 0 )
                {
                    if ( t <= det )  // region 0 (interior)
                    {
                        // minimum at two interior points of 3D lines
                        double invDet = 1.0f/det;
                        s *= invDet;
                        t *= invDet;
                        return DIST(s*(A*s+B*t+2*D)+t*(B*s+C*t+2*E)+F);
                    }
                    else  // region 3 (side)
                    {
                        t = 1;
                        tmp = B+D;
                        if ( tmp >= 0 )
                        {
                            s = 0;
                            return DIST(C+2*E+F);
                        }
                        else if ( -tmp >= A )
                        {
                            s = 1;
                            return DIST(A+C+F+2*(E+tmp));
                        }
                        else
                        {
                            s = -tmp/A;
                            return DIST(tmp*s+C+2*E+F);
                        }
                    }
                }
                else  // region 7 (side)
                {
                    t = 0;
                    if ( D >= 0 )
                    {
                        s = 0;
                        return DIST(F);
                    }
                    else if ( -D >= A )
                    {
                        s = 1;
                        return DIST(A+2*D+F);
                    }
                    else
                    {
                        s = -D/A;
                        return DIST(D*s+F);
                    }
                }
            }
            else
            {
                if ( t >= 0 )
                {
                    if ( t <= det )  // region 1 (side)
                    {
                        s = 1;
                        tmp = B+E;
                        if ( tmp >= 0 )
                        {
                            t = 0;
                            return DIST(A+2*D+F);
                        }
                        else if ( -tmp >= C )
                        {
                            t = 1;
                            return DIST(A+C+F+2*(D+tmp));
                        }
                        else
                        {
                            t = -tmp/C;
                            return DIST(tmp*t+A+2*D+F);
                        }
                    }
                    else  // region 2 (corner)
                    {
                        tmp = B+D;
                        if ( -tmp <= A )
                        {
                            t = 1;
                            if ( tmp >= 0 )
                            {
                                s = 0;
                                return DIST(C+2*E+F);
                            }
                            else
                            {
                                 s = -tmp/A;
                                 return DIST(tmp*s+C+2*E+F);
                            }
                        }
                        else
                        {
                            s = 1;
                            tmp = B+E;
                            if ( tmp >= 0 )
                            {
                                t = 0;
                                return DIST(A+2*D+F);
                            }
                            else if ( -tmp >= C )
                            {
                                t = 1;
                                return DIST(A+C+F+2*(D+tmp));
                            }
                            else
                            {
                                t = -tmp/C;
                                return DIST(tmp*t+A+2*D+F);
                            }
                        }
                    }
                }
                else  // region 8 (corner)
                {
                    if ( -D < A )
                    {
                        t = 0;
                        if ( D >= 0 )
                        {
                            s = 0;
                            return DIST(F);
                        }
                        else
                        {
                            s = -D/A;
                            return DIST(D*s+F);
                        }
                    }
                    else
                    {
                        s = 1;
                        tmp = B+E;
                        if ( tmp >= 0 )
                        {
                            t = 0;
                            return DIST(A+2*D+F);
                        }
                        else if ( -tmp >= C )
                        {
                            t = 1;
                            return DIST(A+C+F+2*(D+tmp));
                        }
                        else
                        {
                            t = -tmp/C;
                            return DIST(tmp*t+A+2*D+F);
                        }
                    }
                }
            }
        }
        else 
        {
            if ( t >= 0 )
            {
                if ( t <= det )  // region 5 (side)
                {
                    s = 0;
                    if ( E >= 0 )
                    {
                        t = 0;
                        return DIST(F);
                    }
                    else if ( -E >= C )
                    {
                        t = 1;
                        return DIST(C+2*E+F);
                    }
                    else
                    {
                        t = -E/C;
                        return DIST(E*t+F);
                    }
                }
                else  // region 4 (corner)
                {
                    tmp = B+D;
                    if ( tmp < 0 )
                    {
                        t = 1;
                        if ( -tmp >= A )
                        {
                            s = 1;
                            return DIST(A+C+F+2*(E+tmp));
                        }
                        else
                        {
                            s = -tmp/A;
                            return DIST(tmp*s+C+2*E+F);
                        }
                    }
                    else
                    {
                        s = 0;
                        if ( E >= 0 )
                        {
                            t = 0;
                            return DIST(F);
                        }
                        else if ( -E >= C )
                        {
                            t = 1;
                            return DIST(C+2*E+F);
                        }
                        else
                        {
                            t = -E/C;
                            return DIST(E*t+F);
                        }
                    }
                }
            }
            else   // region 6 (corner)
            {
                if ( D < 0 )
                {
                    t = 0;
                    if ( -D >= A )
                    {
                        s = 1;
                        return DIST(A+2*D+F);
                    }
                    else
                    {
                        s = -D/A;
                        return DIST(D*s+F);
                    }
                }
                else
                {
                    s = 0;
                    if ( E >= 0 )
                    {
                        t = 0;
                        return DIST(F);
                    }
                    else if ( -E >= C )
                    {
                        t = 1;
                        return DIST(C+2*E+F);
                    }
                    else
                    {
                        t = -E/C;
                        return DIST(E*t+F);
                    }
                }
            }
        }
    }
    else
    {
        // line segments are parallel
        if ( B > 0 )
        {
            // direction vectors form an obtuse angle
            if ( D >= 0 )
            {
                s = 0;
                t = 0;
                return DIST(F);
            }
            else if ( -D <= A )
            {
                s = -D/A;
                t = 0;
                return DIST(D*s+F);
            }
            else
            {
                E = -Dot(seg1.m,diff);
                s = 1;
                tmp = A+D;
                if ( -tmp >= B )
                {
                    t = 1;
                    return DIST(A+C+F+2*(B+D+E));
                }
                else
                {
                    t = -tmp/B;
                    return DIST(A+2*D+F+t*(C*t+2*(B+E)));
                }
            }
        }
        else
        {
            // direction vectors form an acute angle
            if ( -D >= A )
            {
                s = 1;
                t = 0;
                return DIST(A+2*D+F);
            }
            else if ( D <= 0 )
            {
                s = -D/A;
                t = 0;
                return DIST(D*s+F);
            }
            else
            {
                E = -Dot(seg1.m,diff);
                s = 0;
                if ( D >= -B )
                {
                    t = 1;
                    return DIST(C+2*E+F);
                }
                else
                {
                    t = -D/B;
                    return DIST(F+t*(2*E+C*t));
                }
            }
        }
    }
}
double AnalyticGeometryTool::MinLineSegmentTriangle ( const Line3 seg,
const Triangle3 tri,
double &  r,
double &  s,
double &  t 
) [private]

Definition at line 5340 of file AnalyticGeometryTool.cpp.

{
    Point3 diff = tri.b - seg.b;
    double A00 = Dot(seg.m,seg.m);
    double A01 = -Dot(seg.m,tri.e0);
    double A02 = -Dot(seg.m,tri.e1);
    double A11 = Dot(tri.e0,tri.e0);
    double A12 = Dot(tri.e0,tri.e1);
    double A22 = Dot(tri.e1,tri.e1);
    double B0  = -Dot(diff,seg.m);
    double B1  = Dot(diff,tri.e0);
    double B2  = Dot(diff,tri.e1);
    double cof00 = A11*A22-A12*A12;
    double cof01 = A02*A12-A01*A22;
    double cof02 = A01*A12-A02*A11;
    double det = A00*cof00+A01*cof01+A02*cof02;

    Line3 triseg;
    Point3 pt;
    double min, min0, r0, s0, t0;

    if ( Abs(det) >= par_tolerance )
    {
        double cof11 = A00*A22-A02*A02;
        double cof12 = A02*A01-A00*A12;
        double cof22 = A00*A11-A01*A01;
        double invDet = 1.0f/det;
        double rhs0 = -B0*invDet;
        double rhs1 = -B1*invDet;
        double rhs2 = -B2*invDet;

        r = cof00*rhs0+cof01*rhs1+cof02*rhs2;
        s = cof01*rhs0+cof11*rhs1+cof12*rhs2;
        t = cof02*rhs0+cof12*rhs1+cof22*rhs2;

        if ( r < 0 )
        {
            if ( s+t <= 1 )
            {
                if ( s < 0 )
                {
                    if ( t < 0 )  // region 4m
                    {
                        // min on face s=0 or t=0 or r=0
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                        triseg.b = tri.b;
                        triseg.m = tri.e0;
                        min0 = MinLineSegmentLineSegment(seg,triseg,r0,s0);
                        t0 = 0;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                        min0 = MinPointTriangle(seg.b,tri,s0,t0);
                        r0 = 0;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                    }
                    else  // region 3m
                    {
                        // min on face s=0 or r=0
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                        min0 = MinPointTriangle(seg.b,tri,s0,t0);
                        r0 = 0;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                    }
                }
                else if ( t < 0 )  // region 5m
                {
                    // min on face t=0 or r=0
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                    min0 = MinPointTriangle(seg.b,tri,s0,t0);
                    r0 = 0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else  // region 0m
                {
                    // min face on r=0
                    min = MinPointTriangle(seg.b,tri,s,t);
                    r = 0;
                }
            }
            else
            {
                if ( s < 0 )  // region 2m
                {
                    // min on face s=0 or s+t=1 or r=0
                    triseg.b = tri.b;
                    triseg.m = tri.e1;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                    min0 = MinPointTriangle(seg.b,tri,s0,t0);
                    r0 = 0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else if ( t < 0 )  // region 6m
                {
                    // min on face t=0 or s+t=1 or r=0
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                    min0 = MinPointTriangle(seg.b,tri,s0,t0);
                    r0 = 0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else  // region 1m
                {
                    // min on face s+t=1 or r=0
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 1-t;
                    min0 = MinPointTriangle(seg.b,tri,s0,t0);
                    r0 = 0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
            }
        }
        else if ( r <= 1 )
        {
            if ( s+t <= 1 )
            {
                if ( s < 0 )
                {
                    if ( t < 0 )  // region 4
                    {
                        // min on face s=0 or t=0
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                        triseg.b = tri.b;
                        triseg.m = tri.e0;
                        min0 = MinLineSegmentLineSegment(seg,triseg,r0,s0);
                        t0 = 0;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                    }
                    else  // region 3
                    {
                        // min on face s=0
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                    }
                }
                else if ( t < 0 )  // region 5
                {
                    // min on face t=0
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                }
                else  // region 0
                {
                    // global minimum is interior, done
                    min = Sqrt(Abs(r*(A00*r+A01*s+A02*t+2.0f*B0)
                          +s*(A01*r+A11*s+A12*t+2.0f*B1)
                          +t*(A02*r+A12*s+A22*t+2.0f*B2)
                          +Dot(diff,diff)));
                }
            }
            else
            {
                if ( s < 0 )  // region 2
                {
                    // min on face s=0 or s+t=1
                    triseg.b = tri.b;
                    triseg.m = tri.e1;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else if ( t < 0 )  // region 6
                {
                    // min on face t=0 or s+t=1
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else  // region 1
                {
                    // min on face s+t=1
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 1-t;
                }
            }
        }
        else  // r > 1
        {
            if ( s+t <= 1 )
            {
                if ( s < 0 )
                {
                    if ( t < 0 )  // region 4p
                    {
                        // min on face s=0 or t=0 or r=1
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                        triseg.b = tri.b;
                        triseg.m = tri.e0;
                        min0 = MinLineSegmentLineSegment(seg,triseg,r0,s0);
                        t0 = 0;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                        pt = seg.b+seg.m;
                        min0 = MinPointTriangle(pt,tri,s0,t0);
                        r0 = 1;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                    }
                    else  // region 3p
                    {
                        // min on face s=0 or r=1
                        triseg.b = tri.b;
                        triseg.m = tri.e1;
                        min = MinLineSegmentLineSegment(seg,triseg,r,t);
                        s = 0;
                        pt = seg.b+seg.m;
                        min0 = MinPointTriangle(pt,tri,s0,t0);
                        r0 = 1;
                        if ( min0 < min )
                        {
                            min = min0;
                            r = r0;
                            s = s0;
                            t = t0;
                        }
                    }
                }
                else if ( t < 0 )  // region 5p
                {
                    // min on face t=0 or r=1
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                    pt = seg.b+seg.m;
                    min0 = MinPointTriangle(pt,tri,s0,t0);
                    r0 = 1;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else  // region 0p
                {
                    // min face on r=1
                    pt = seg.b+seg.m;
                    min = MinPointTriangle(pt,tri,s,t);
                    r = 1;
                }
            }
            else
            {
                if ( s < 0 )  // region 2p
                {
                    // min on face s=0 or s+t=1 or r=1
                    triseg.b = tri.b;
                    triseg.m = tri.e1;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                    pt = seg.b+seg.m;
                    min0 = MinPointTriangle(pt,tri,s0,t0);
                    r0 = 1;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else if ( t < 0 )  // region 6p
                {
                    // min on face t=0 or s+t=1 or r=1
                    triseg.b = tri.b;
                    triseg.m = tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,s);
                    t = 0;
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
                    s0 = 1-t0;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                    pt = seg.b+seg.m;
                    min0 = MinPointTriangle(pt,tri,s0,t0);
                    r0 = 1;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
                else  // region 1p
                {
                    // min on face s+t=1 or r=1
                    triseg.b = tri.b+tri.e0;
                    triseg.m = tri.e1-tri.e0;
                    min = MinLineSegmentLineSegment(seg,triseg,r,t);
                    s = 1-t;
                    pt = seg.b+seg.m;
                    min0 = MinPointTriangle(pt,tri,s0,t0);
                    r0 = 1;
                    if ( min0 < min )
                    {
                        min = min0;
                        r = r0;
                        s = s0;
                        t = t0;
                    }
                }
            }
        }
    }
    else
    {
        // line and triangle are parallel
        triseg.b = tri.b;
        triseg.m = tri.e0;
        min = MinLineSegmentLineSegment(seg,triseg,r,s);
        t = 0;

        triseg.m = tri.e1;
        min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
        s0 = 0;
        if ( min0 < min )
        {
            min = min0;
            r = r0;
            s = s0;
            t = t0;
        }

        triseg.b = triseg.b + tri.e0;
        triseg.m = triseg.m - tri.e0;
        min0 = MinLineSegmentLineSegment(seg,triseg,r0,t0);
        s0 = 1-t0;
        if ( min0 < min )
        {
            min = min0;
            r = r0;
            s = s0;
            t = t0;
        }

        min0 = MinPointTriangle(seg.b,tri,s0,t0);
        r0 = 0;
        if ( min0 < min )
        {
            min = min0;
            r = r0;
            s = s0;
            t = t0;
        }

        pt = seg.b+seg.m;
        min0 = MinPointTriangle(pt,tri,s0,t0);
        r0 = 1;
        if ( min0 < min )
        {
            min = min0;
            r = r0;
            s = s0;
            t = t0;
        }
    }

    return min;
}
double AnalyticGeometryTool::MinPointTriangle ( const Point3 p,
const Triangle3 tri,
double &  s,
double &  t 
)

Finds minimum distance beween a triange and point in 3D (MAGIC)

Definition at line 5015 of file AnalyticGeometryTool.cpp.

{
    Point3 diff = tri.b - p;
    double A = Dot(tri.e0,tri.e0);
    double B = Dot(tri.e0,tri.e1);
    double C = Dot(tri.e1,tri.e1);
    double D = Dot(tri.e0,diff);
    double E = Dot(tri.e1,diff);
    double F = Dot(diff,diff);
    double det = Abs(A*C-B*B);  // A*C-B*B = |Cross(e0,e1)|^2 >= 0
    // A-2*B+C = Dot(e0,e0)-2*Dot(e0,e1)+Dot(e1,e1) = |e0-e1|^2 > 0

    s = B*E-C*D;
    t = B*D-A*E;

    if ( s+t <= det )
    {
        if ( s < 0 )
        {
            if ( t < 0 )  // region 4
            {
                if ( D < 0 )
                {
                    t = 0;
                    if ( -D >= A )
                    {
                        s = 1;
                        return DIST(A+2*D+F);
                    }
                    else
                    {
                        s = -D/A;
                        return DIST(D*s+F);
                    }
                }
                else
                {
                    s = 0;
                    if ( E >= 0 )
                    {
                        t = 0;
                        return DIST(F);
                    }
                    else if ( -E >= C )
                    {
                        t = 1;
                        return DIST(C+2*E+F);
                    }
                    else
                    {
                        t = -E/C;
                        return DIST(E*t+F);
                    }
                }
            }
            else  // region 3
            {
                s = 0;
                if ( E >= 0 )
                {
                    t = 0;
                    return DIST(F);
                }
                else if ( -E >= C )
                {
                    t = 1;
                    return DIST(C+2*E+F);
                }
                else
                {
                    t = -E/C;
                    return DIST(E*t+F);
                }
            }
        }
        else if ( t < 0 )  // region 5
        {
            t = 0;
            if ( D >= 0 )
            {
                s = 0;
                return DIST(F);
            }
            else if ( -D >= A )
            {
                s = 1;
                return DIST(A+2*D+F);
            }
            else
            {
                s = -D/A;
                return DIST(D*s+F);
            }
        }
        else  // region 0
        {
            // minimum at interior point
          if( det == 0.0 )
          {
            //PRINT_WARNING( "Found zero determinant\n" );
            return CUBIT_DBL_MAX;
          }

            double invDet = 1.0f/det;
            s *= invDet;
            t *= invDet;
            return DIST(s*(A*s+B*t+2*D)+t*(B*s+C*t+2*E)+F);
        }
    }
    else
    {
        double tmp0, tmp1, numer, denom;

        if ( s < 0 )  // region 2
        {
            tmp0 = B+D;
            tmp1 = C+E;
            if ( tmp1 > tmp0 )
            {
                numer = tmp1 - tmp0;
                denom = A-2*B+C;
                if ( numer >= denom )
                {
                    s = 1;
                    t = 0;
                    return DIST(A+2*D+F);
                }
                else
                {
                    s = numer/denom;
                    t = 1-s;
                    return DIST(s*(A*s+B*t+2*D)+t*(B*s+C*t+2*E)+F);
                }
            }
            else
            {
                s = 0;
                if ( tmp1 <= 0 )
                {
                    t = 1;
                    return DIST(C+2*E+F);
                }
                else if ( E >= 0 )
                {
                    t = 0;
                    return DIST(F);
                }
                else
                {
                    t = -E/C;
                    return DIST(E*t+F);
                }
            }
        }
        else if ( t < 0 )  // region 6
        {
            tmp0 = B+E;
            tmp1 = A+D;
            if ( tmp1 > tmp0 )
            {
                numer = tmp1 - tmp0;
                denom = A-2*B+C;
                if ( numer >= denom )
                {
                    t = 1;
                    s = 0;
                    return DIST(C+2*E+F);
                }
                else
                {
                    t = numer/denom;
                    s = 1-t;
                    return DIST(s*(A*s+B*t+2*D)+t*(B*s+C*t+2*E)+F);
                }
            }
            else
            {
                t = 0;
                if ( tmp1 <= 0 )
                {
                    s = 1;
                    return DIST(A+2*D+F);
                }
                else if ( D >= 0 )
                {
                    s = 0;
                    return DIST(F);
                }
                else
                {
                    s = -D/A;
                    return DIST(D*s+F);
                }
            }
        }
        else  // region 1
        {
            numer = C+E-B-D;
            if ( numer <= 0 )
            {
                s = 0;
                t = 1;
                return DIST(C+2*E+F);
            }
            else
            {
                denom = A-2*B+C;
                if ( numer >= denom )
                {
                    s = 1;
                    t = 0;
                    return DIST(A+2*D+F);
                }
                else
                {
                    s = numer/denom;
                    t = 1-s;
                    return DIST(s*(A*s+B*t+2*D)+t*(B*s+C*t+2*E)+F);
                }
            }
        }
    }
}
double AnalyticGeometryTool::MinTriangleTriangle ( Triangle3 tri0,
Triangle3 tri1,
double &  s,
double &  t,
double &  u,
double &  v 
)

Finds minimum distance between two triangles in 3D (MAGIC)

Definition at line 4412 of file AnalyticGeometryTool.cpp.

{
    double s0, t0, u0, v0, min, min0;
    Line3 line;

    // compare edges of tri0 against all of tri1
    line.b = tri0.b;
    line.m = tri0.e0;
    min = MinLineSegmentTriangle(line,tri1,s,u,v);
    t = 0;

    line.m = tri0.e1;
    min0 = MinLineSegmentTriangle(line,tri1,t0,u0,v0);
    s0 = 0;
    if ( min0 < min )
    {
        min = min0;
        s = s0;
        t = t0;
        u = u0;
        v = v0;
    }

    line.b = line.b + tri0.e0;
    line.m = line.m - tri0.e0;
    min0 = MinLineSegmentTriangle(line,tri1,t0,u0,v0);
    s0 = 1-t0;
    if ( min0 < min )
    {
        min = min0;
        s = s0;
        t = t0;
        u = u0;
        v = v0;
    }

    // compare edges of tri1 against all of tri0
    line.b = tri1.b;
    line.m = tri1.e0;
    min0 = MinLineSegmentTriangle(line,tri0,u0,s0,t0);
    v0 = 0;
    if ( min0 < min )
    {
        min = min0;
        s = s0;
        t = t0;
        u = u0;
        v = v0;
    }

    line.m = tri1.e1;
    min0 = MinLineSegmentTriangle(line,tri0,v0,s0,t0);
    u0 = 0;
    if ( min0 < min )
    {
        min = min0;
        s = s0;
        t = t0;
        u = u0;
        v = v0;
    }

    line.b = line.b + tri1.e0;
    line.m = line.m - tri1.e0;
    min0 = MinLineSegmentTriangle(line,tri0,v0,s0,t0);
    u0 = 1-v0;
    if ( min0 < min )
    {
        min = min0;
        s = s0;
        t = t0;
        u = u0;
        v = v0;
    }

    return min;
}
CubitStatus AnalyticGeometryTool::mirror_pnt ( double  pnt[3],
double  pln_orig[3],
double  pln_norm[3],
double  m_pnt[3] 
)

Function to mirror a point about a plane. The mirror point is the same distance from the plane as the original, but on the opposite side of the plane. Function returns CUBIT_FAILURE if the point is on the plane - in which case the point is just copied.

Definition at line 1040 of file AnalyticGeometryTool.cpp.

{
   double int_pnt[3],
      vec[3];
   
   // Find intersection of point and plane 
   if (int_pnt_pln(pnt, pln_orig, pln_norm, int_pnt)) {
      // If intersection is on the plane, return 
      copy_pnt(pnt, m_pnt);
      return CUBIT_FAILURE;
   }
   
   // Find vector from pnt to int_pnt 
   get_vec(pnt, int_pnt, vec);
   
   // Traverse twice the length of vec in vec direction 
   m_pnt[0] = pnt[0] + 2.0 * vec[0];
   m_pnt[1] = pnt[1] + 2.0 * vec[1];
   m_pnt[2] = pnt[2] + 2.0 * vec[2];
   
   return CUBIT_SUCCESS;   
}
void AnalyticGeometryTool::mtx_to_angs ( double  mtx3x3[3][3],
double &  ax,
double &  ay,
double &  az 
)

Gets rotation angles to rotate one system to another - returned rotation angles are about global system. To use the result from this function, align the unoriented object's origin with the global origin, then apply the rotation angles returned from this routine to the object in the order of ax,ay,az about the global origin. The object will be oriented such that its xyz axes will point in the same direction as the object whose transformation matrix was given. The object can then be translated.

Definition at line 446 of file AnalyticGeometryTool.cpp.

{   
//   METHOD:
//   o Rotate x-vector onto xz plane
//   o  Check xp dotted into y  
//   o   If xp dot y is zero  ==> az = 0 (x-vector already in xz plane)
//   o   Otherwise, compute rotation of vector into xz plane to acquire *az     
//   o    Use atan2 (on x-vector) to get *az
//   o    Rotate the system about z 
//   o Use atan2 function (on x-vector in xz plane) to determine *ay 
//   o Rotate the system about y      
//   o Compute ax using y-vector
//   o Resultant angles are negated (to reverse above procedure)
 
   double x[3];                    // x-axis vector 
   double y[3];                    // y-axis vector 
   double z[3];                    // z-axis vector 
   double ar[3][3];                // Rotation matrix 
   double sinr,cosr;               // Used for atan2 function 
   double work_sys[3][3];          // Temporary holder for system 
   double *xp = work_sys[0];       // x-axis vector of system: x-primed 
   double *yp = work_sys[1];       // y-axis vector of system: y-primed 
   
   x[0] = 1.0; x[1] = 0.0; x[2] = 0.0;
   y[0] = 0.0; y[1] = 1.0; y[2] = 0.0;
   z[0] = 0.0; z[1] = 0.0; z[2] = 1.0;
   
   if (!mtx3x3)
      return;
   
   // Copy matrix over to work csys 
   copy_mtx(mtx3x3,work_sys);
      
   // Check xp dotted into y 
   //   If xp dot y is zero  ==> az = 0 
   
   // Otherwise, compute rotation of vector into xz plane to acquire *az 
   
   if (dbl_eq(dot_vec(xp,y), 0.0))
      
      az = 0.0;
     
   else {
      
      /* 
      Compute *az - rotate xp to xz-plane about z-axis          
               y    xp                                           
               |  /                                              
               | /  negative angle about z (negative of atan2)   
                -----x              (use RH rule about z)        
                 \                                               
                  \ positive angle about z (negative of atan2)         
                   xp                                            
      */
     
      sinr = dot_vec(xp,y);
      cosr = dot_vec(xp,x);
      az = - atan2(sinr, cosr);
      
      // Rotate the system about z 
      create_rotation_mtx(az,z,ar);
      rotate_system(ar,work_sys,work_sys);
      
   }
   
      /* 
      Compute *ay - rotate xp to x-axis about y-axis            
              z    xp                                           
              |  /                                              
              | /  positive angle about y (positive of atan2)   
               -----x     (use RH rule about y)                 
                \                                               
                 \ negative angle about y (positive of atan2)   
                  xp                                            
      */
   
   sinr = dot_vec(xp,z);
   cosr = dot_vec(xp,x);
   ay = atan2(sinr, cosr);
   
   // Rotate the system about y       
   create_rotation_mtx(ay,y,ar);      
   rotate_system(ar,work_sys,work_sys);
  
   /*
      Compute *ax - rotate yp to y-axis about x-axis            
              z    yp                                           
              |  /                                              
              | /  negative angle about x (negative of atan2)   
               -----y     (use RH rule about x)                 
                \                                               
                 \ positive angle about x (negative of atan2)   
                  yp                                            
   */
   
   sinr = dot_vec(yp,z);
   cosr = dot_vec(yp,y);
   ax = atan2(sinr,cosr);     // Negative of negative - see below 
   
   // Negate above angles for rotation of the system back to original 
   az = -(az);
   ay = -(ay);
   
   // Make sure near zero angles are actually zero 
   if (dbl_eq(ax, 0.0))
      ax = 0.0;
   
   if (dbl_eq(ay, 0.0))
      ay = 0.0;   
}
void AnalyticGeometryTool::mtx_to_angs ( double  mtx4x4[4][4],
double &  ax,
double &  ay,
double &  az 
)

Gets rotation angles to rotate one system to another - returned rotation angles are about global system. To use the result from this function, align the unoriented object's origin with the global origin, then apply the rotation angles returned from this routine to the object in the order of ax,ay,az about the global origin. The object will be oriented such that its xyz axes will point in the same direction as the object whose transformation matrix was given. The object can then be translated.

Definition at line 559 of file AnalyticGeometryTool.cpp.

{   
   double work_sys[3][3];

   if(!mtx4x4)
      return;

   copy_mtx(mtx4x4,work_sys);
   mtx_to_angs( work_sys, ax, ay, az );
}
void AnalyticGeometryTool::mult_mtx ( double  a[3][3],
double  b[3][3],
double  d[3][3] 
)

Definition at line 707 of file AnalyticGeometryTool.cpp.

{       
   double c[3][3];    // working buffer 
      
   if( a != 0 && b != 0 ) {   // a & b are valid 
            
      c[0][0] = (  a[0][0] * b[0][0] + a[0][1] * b[1][0] 
                 + a[0][2] * b[2][0]);
      c[1][0] = (  a[1][0] * b[0][0] + a[1][1] * b[1][0] 
                 + a[1][2] * b[2][0]);
      c[2][0] = (  a[2][0] * b[0][0] + a[2][1] * b[1][0] 
                 + a[2][2] * b[2][0]);
      
      c[0][1] = (  a[0][0] * b[0][1] + a[0][1] * b[1][1] 
                 + a[0][2] * b[2][1]);
      c[1][1] = (  a[1][0] * b[0][1] + a[1][1] * b[1][1] 
                 + a[1][2] * b[2][1]);
      c[2][1] = (  a[2][0] * b[0][1] + a[2][1] * b[1][1] 
                 + a[2][2] * b[2][1]);
      
      c[0][2] = (  a[0][0] * b[0][2] + a[0][1] * b[1][2] 
                 + a[0][2] * b[2][2]);
      c[1][2] = (  a[1][0] * b[0][2] + a[1][1] * b[1][2]
                 + a[1][2] * b[2][2]);
      c[2][2] = (  a[2][0] * b[0][2] + a[2][1] * b[1][2]
                 + a[2][2] * b[2][2]);
      
      copy_mtx(c, d);
      
   }
   else if (a) {                 // b equals 0     
      copy_mtx(a, d);    
   }
   else if (b) {                  // a equals 0 
      copy_mtx(b, d);    
   }
   else {                         // a & b equal 0 
      
      copy_mtx(AGT_IDENTITY_MTX_3X3, d);     
   }
}
void AnalyticGeometryTool::mult_mtx ( double  a[4][4],
double  b[4][4],
double  d[4][4] 
)

Definition at line 750 of file AnalyticGeometryTool.cpp.

{       
   double c[4][4];    // working buffer 
      
   if( a != 0 && b != 0 ) {   // a & b are valid 
            
      c[0][0] = (  a[0][0] * b[0][0] + a[0][1] * b[1][0] 
                 + a[0][2] * b[2][0] + a[0][3] * b[3][0]);
      c[1][0] = (  a[1][0] * b[0][0] + a[1][1] * b[1][0] 
                 + a[1][2] * b[2][0] + a[1][3] * b[3][0]);
      c[2][0] = (  a[2][0] * b[0][0] + a[2][1] * b[1][0] 
                 + a[2][2] * b[2][0] + a[2][3] * b[3][0]);
      c[3][0] = (  a[3][0] * b[0][0] + a[3][1] * b[1][0] 
                 + a[3][2] * b[2][0] + a[3][3] * b[3][0]);
      
      c[0][1] = (  a[0][0] * b[0][1] + a[0][1] * b[1][1] 
                 + a[0][2] * b[2][1] + a[0][3] * b[3][1]);
      c[1][1] = (  a[1][0] * b[0][1] + a[1][1] * b[1][1] 
                 + a[1][2] * b[2][1] + a[1][3] * b[3][1]);
      c[2][1] = (  a[2][0] * b[0][1] + a[2][1] * b[1][1] 
                 + a[2][2] * b[2][1] + a[2][3] * b[3][1]);
      c[3][1] = (  a[3][0] * b[0][1] + a[3][1] * b[1][1] 
                 + a[3][2] * b[2][1] + a[3][3] * b[3][1]);
      
      c[0][2] = (  a[0][0] * b[0][2] + a[0][1] * b[1][2] 
                 + a[0][2] * b[2][2] + a[0][3] * b[3][2]);
      c[1][2] = (  a[1][0] * b[0][2] + a[1][1] * b[1][2]
                 + a[1][2] * b[2][2] + a[1][3] * b[3][2]);
      c[2][2] = (  a[2][0] * b[0][2] + a[2][1] * b[1][2]
                 + a[2][2] * b[2][2] + a[2][3] * b[3][2]);
      c[3][2] = (  a[3][0] * b[0][2] + a[3][1] * b[1][2]
                 + a[3][2] * b[2][2] + a[3][3] * b[3][2]);
      
      c[0][3] = (  a[0][0] * b[0][3] + a[0][1] * b[1][3]
                 + a[0][2] * b[2][3] + a[0][3] * b[3][3]);
      c[1][3] = (  a[1][0] * b[0][3] + a[1][1] * b[1][3]
                 + a[1][2] * b[2][3] + a[1][3] * b[3][3]);
      c[2][3] = (  a[2][0] * b[0][3] + a[2][1] * b[1][3]
                 + a[2][2] * b[2][3] + a[2][3] * b[3][3]);
      c[3][3] = (  a[3][0] * b[0][3] + a[3][1] * b[1][3]
                 + a[3][2] * b[2][3] + a[3][3] * b[3][3]);
      
      copy_mtx(c, d);        
   }
   else if (a) {                 // b equals 0     
      copy_mtx(a, d);    
   }
   else if (b) {                  // a equals 0 
      copy_mtx(b, d);    
   }
   else {                         // a & b equal 0 
      
      copy_mtx(AGT_IDENTITY_MTX_4X4, d);     
   }
}
void AnalyticGeometryTool::mult_vecxconst ( double  constant,
double  vec[3],
double  vec_out[3] 
)

Multiples a vector by a constant (vec_out can equal vec).

Definition at line 1237 of file AnalyticGeometryTool.cpp.

{
   // multiply each element of the vector by the constant
   vec_out[0] = constant * vec[0];
   vec_out[1] = constant * vec[1];
   vec_out[2] = constant * vec[2];   
}
CubitStatus AnalyticGeometryTool::next_pnt ( double  str_pnt[3],
double  vec_dir[3],
double  len,
double  new_pnt[3] 
)

Given start pnt, vec dir and length find next point in space. The vector direction is first unitized then the following formula is used: new_point[] = start_point[] + (length * unit_vector[]) new_pnt can be the same as input point (overwrites old location).

Definition at line 1067 of file AnalyticGeometryTool.cpp.

{        
   double uv[3];  // unit vector representation of vector direction 
   
   // unitize specified direction 
   if (!unit_vec(vec_dir,uv)) {
      copy_pnt(str_pnt, new_pnt);
      return CUBIT_FAILURE;
   }
   
   // determine next point in space 
   
   new_pnt[0] = str_pnt[0] + (len * uv[0]);     
   new_pnt[1] = str_pnt[1] + (len * uv[1]);     
   new_pnt[2] = str_pnt[2] + (len * uv[2]); 
   
   return CUBIT_SUCCESS;
}
void AnalyticGeometryTool::Normal ( Triangle3 tri,
double  normal[3] 
)

Finds normal vector of given triangle.

Definition at line 6217 of file AnalyticGeometryTool.cpp.

{
  double vec1[3];
  vec1[0]=tri.e0.x; vec1[1]=tri.e0.y; vec1[2]=tri.e0.z;

  double vec2[3];
  vec2[0]=tri.e1.x; vec2[1]=tri.e1.y; vec2[2]=tri.e1.z;

  cross_vec_unit( vec1, vec2, normal );
}
void AnalyticGeometryTool::orth_vecs ( double  uvect[3],
double  vvect[3],
double  wvect[3] 
)

Finds 2 arbitrary orthoganal vectors to the first.

Definition at line 1148 of file AnalyticGeometryTool.cpp.

{
   double x[3];
   unsigned short i = 0;
   unsigned short imin = 3;
   double rmin = 1.0E20;
   unsigned short iperm1[3];
   unsigned short iperm2[3];
   unsigned short cont_flag = 1;
   double vec[3];
   
   // Initialize perm flags
   iperm1[0] = 1; iperm1[1] = 2; iperm1[2] = 0;
   iperm2[0] = 2; iperm2[1] = 0; iperm2[2] = 1;
   
   // unitize vector 
   
   unit_vec(uvect,vec);
   
   while (i<3 && cont_flag ) {
      if (dbl_eq(vec[i], 0.0)) {
         vvect[i] = 1.0;
         vvect[iperm1[i]] = 0.0;
         vvect[iperm2[i]] = 0.0;
         cont_flag = 0;
      }
      
      if (fabs(vec[i]) < rmin) {
         imin = i;
         rmin = fabs(vec[i]);
      }
      ++i;
   }
   
   if (cont_flag) {
      x[imin] = 1.0;
      x[iperm1[imin]] = 0.0;
      x[iperm2[imin]] = 0.0;
      
      // determine cross product 
      cross_vec_unit(vec,x,vvect);
   }
   
   // cross vector to determine last orthogonal vector 
   cross_vec(vvect,vec,wvect);
}
CubitBoolean AnalyticGeometryTool::pnt_eq ( double  pnt1[3],
double  pnt2[3] 
)

Compares two points to determine if they are equivalent. The equivalence is based on epsilon (get_epsilon, set_epsilon). If the distance between them is less than epsilon, the points are considered equal.

Definition at line 1030 of file AnalyticGeometryTool.cpp.

{
   double x = pnt2[0] - pnt1[0];  // difference in the x direction 
   double y = pnt2[1] - pnt1[1];  // difference in the y direction 
   double z = pnt2[2] - pnt1[2];  // difference in the z direction 
  
   return (dbl_eq(sqrt(x*x + y*y + z*z), 0.0));
}
void AnalyticGeometryTool::print_mtx ( double  mtx[3][3])

Prints matrix values, for debugging.

Definition at line 980 of file AnalyticGeometryTool.cpp.

{
   PRINT_INFO( "%f %f %f\n", mtx[0][0], mtx[0][1], mtx[0][2] );
   PRINT_INFO( "%f %f %f\n", mtx[1][0], mtx[1][1], mtx[1][2] );
   PRINT_INFO( "%f %f %f\n", mtx[2][0], mtx[2][1], mtx[2][2] );
}
void AnalyticGeometryTool::print_mtx ( double  mtx[4][4])

Prints matrix values, for debugging.

Definition at line 987 of file AnalyticGeometryTool.cpp.

{
   PRINT_INFO( "%f %f %f %f\n", mtx[0][0], mtx[0][1], mtx[0][2], mtx[0][3] );
   PRINT_INFO( "%f %f %f %f\n", mtx[1][0], mtx[1][1], mtx[1][2], mtx[1][3] );
   PRINT_INFO( "%f %f %f %f\n", mtx[2][0], mtx[2][1], mtx[2][2], mtx[2][3] );
   PRINT_INFO( "%f %f %f %f\n", mtx[3][0], mtx[3][1], mtx[3][2], mtx[3][3] );
}
double AnalyticGeometryTool::ProjectedOverlap ( Triangle3 tri1,
Triangle3 tri2,
bool  draw_overlap = false 
)

Projects tri2 to the plane of tri1 and finds the overlap area.

Definition at line 6229 of file AnalyticGeometryTool.cpp.

{
  // Transform both triangles into local coordinate system of tri1
  // to eliminate z-coordinate. 

  // Use normal as z-axis
  double z[3];
  Normal( tri1, z );

  // x-axis goes from b to e0
  double x[3];
  x[0] = tri1.e0.x;
  x[1] = tri1.e0.y;
  x[2] = tri1.e0.z;

  // Cross z with x to get y-axis
  double y[3];
  cross_vec( z, x, y );

  // Unitize them all
  unit_vec( x, x );
  unit_vec( y, y );
  unit_vec( z, z );

  double origin[3];
  origin[0] = tri1.b.x;
  origin[1] = tri1.b.y;
  origin[2] = tri1.b.z;

  double mtxTriOne2Global[4][4];
  vecs_to_mtx( x, y, z, origin, mtxTriOne2Global ); // Really mtxTriOne2Global
  double mtxGlobal2TriOne[4][4];
  inv_trans_mtx( mtxTriOne2Global, mtxGlobal2TriOne );
  
  double v0[3], v1[3], v2[3]; // Vertices of triangle
  v0[0]=tri1.b.x;   v0[1]=tri1.b.y;  v0[2]=tri1.b.z;
  v1[0]=tri1.e0.x+tri1.b.x; v1[1]=tri1.e0.y+tri1.b.y; v1[2]=tri1.e0.z+tri1.b.z;
  v2[0]=tri1.e1.x+tri1.b.x; v2[1]=tri1.e1.y+tri1.b.y; v2[2]=tri1.e1.z+tri1.b.z;
  transform_pnt( mtxGlobal2TriOne, v0, v0 );
  transform_pnt( mtxGlobal2TriOne, v1, v1 );
  transform_pnt( mtxGlobal2TriOne, v2, v2 );

  // Load this into a 2D triangle T1 ( z-axis is now zero )
  Triangle T1;
  T1.v[0].x = v0[0]; T1.v[0].y = v0[1];
  T1.v[1].x = v1[0]; T1.v[1].y = v1[1]; 
  T1.v[2].x = v2[0]; T1.v[2].y = v2[1];

  // Setup the next triangle, in this coordinate system
  v0[0]=tri2.b.x;   v0[1]=tri2.b.y;  v0[2]=tri2.b.z;
  v1[0]=tri2.e0.x+tri2.b.x; v1[1]=tri2.e0.y+tri2.b.y; v1[2]=tri2.e0.z+tri2.b.z;
  v2[0]=tri2.e1.x+tri2.b.x; v2[1]=tri2.e1.y+tri2.b.y; v2[2]=tri2.e1.z+tri2.b.z;

  // The following creates coordinates of tri2 in tri1's csys
  transform_pnt( mtxGlobal2TriOne, v0, v0 );
  transform_pnt( mtxGlobal2TriOne, v1, v1 );
  transform_pnt( mtxGlobal2TriOne, v2, v2 );

  // Now we can project tri2 to tri1 simply by dropping the z-coordinate
  Triangle T2;
  T2.v[0].x = v0[0]; T2.v[0].y = v0[1];
  T2.v[1].x = v1[0]; T2.v[1].y = v1[1]; 
  T2.v[2].x = v2[0]; T2.v[2].y = v2[1];

  // Now find area of overlap
  if( draw_overlap )
  {
    AgtTriList* list = Intersection(&T1, &T2);
    double overlap_area = Area(list);
    while ( list )
    {
      Triangle *tmp_tri = list->tri;
      
      v0[0]=tmp_tri->v[0].x;   
      v0[1]=tmp_tri->v[0].y;   
      v0[2]=0;

      v1[0]=tmp_tri->v[1].x;
      v1[1]=tmp_tri->v[1].y;
      v1[2]=0;
      
      v2[0]=tmp_tri->v[2].x; 
      v2[1]=tmp_tri->v[2].y; 
      v2[2]=0;
      
      transform_pnt( mtxTriOne2Global, v0, v0 );
      transform_pnt( mtxTriOne2Global, v1, v1 );
      transform_pnt( mtxTriOne2Global, v2, v2 );

      //Now Draw the thing
       GMem g_mem;
       g_mem.allocate_tri(1);
       g_mem.pointListCount = 3;
      
       g_mem.point_list()[0].x = v0[0]; 
       g_mem.point_list()[0].y = v0[1]; 
       g_mem.point_list()[0].z = v0[2]; 
        
       g_mem.point_list()[1].x = v1[0]; 
       g_mem.point_list()[1].y = v1[1]; 
       g_mem.point_list()[1].z = v1[2]; 
        
       g_mem.point_list()[2].x = v2[0]; 
       g_mem.point_list()[2].y = v2[1]; 
       g_mem.point_list()[2].z = v2[2]; 

       g_mem.facet_list()[0] = 3;
       g_mem.facet_list()[1] = 0;
       g_mem.facet_list()[2] = 1;
       g_mem.facet_list()[3] = 2;

       g_mem.fListCount = 1;
       GfxPreview::draw_polygon(g_mem.point_list(),3,4,4,true);

      AgtTriList* save = list;
      list = list->next;
      delete save->tri;
      delete save;
    }
    delete list;
    return overlap_area;
  }
  else
    return AreaIntersection( T1, T2 );
}
void AnalyticGeometryTool::reverse_vec ( double  vin[3],
double  vout[3] 
)

Multiples -1.0 by a vector (vout can equal vin).

Definition at line 1248 of file AnalyticGeometryTool.cpp.

{
   // Multiply the vector components by a -1.0    
   mult_vecxconst(-1.0, vin, vout);      
}
void AnalyticGeometryTool::rotate_system ( double  mtx[3][3],
double  sys_in[3][3],
double  sys_out[3][3] 
)

Definition at line 572 of file AnalyticGeometryTool.cpp.

{
   double sys_tmp[3][3];
   double *p_sys_tmp;
   
   // Check to see if rotating in place 
   if (sys_in == sys_out) {
      copy_mtx( sys_in, sys_tmp );
      p_sys_tmp = (double *)sys_tmp; 
   }
   else 
      // Have sys_tmp point at outgoing memory location 
      p_sys_tmp = (double *)sys_out;
   
   
   // X-vector 
   p_sys_tmp[0] =   mtx[0][0] * sys_in[0][0]
                  + mtx[1][0] * sys_in[0][1]
                  + mtx[2][0] * sys_in[0][2];
   p_sys_tmp[1] =   mtx[0][1] * sys_in[0][0]
                  + mtx[1][1] * sys_in[0][1]
                  + mtx[2][1] * sys_in[0][2];
   p_sys_tmp[2] =   mtx[0][2] * sys_in[0][0]
                  + mtx[1][2] * sys_in[0][1]
                  + mtx[2][2] * sys_in[0][2];   
   
   // Y-vector 
   p_sys_tmp[3] =   mtx[0][0] * sys_in[1][0]
                  + mtx[1][0] * sys_in[1][1]
                  + mtx[2][0] * sys_in[1][2];
   p_sys_tmp[4] =   mtx[0][1] * sys_in[1][0]
                  + mtx[1][1] * sys_in[1][1]
                  + mtx[2][1] * sys_in[1][2];
   p_sys_tmp[5] =   mtx[0][2] * sys_in[1][0]
                  + mtx[1][2] * sys_in[1][1]
                  + mtx[2][2] * sys_in[1][2];
   
   // Z-vector 
   p_sys_tmp[6] =   mtx[0][0] * sys_in[2][0]
                  + mtx[1][0] * sys_in[2][1]
                  + mtx[2][0] * sys_in[2][2];
   p_sys_tmp[7] =   mtx[0][1] * sys_in[2][0]
                  + mtx[1][1] * sys_in[2][1]
                  + mtx[2][1] * sys_in[2][2];
   p_sys_tmp[8] =   mtx[0][2] * sys_in[2][0]
                  + mtx[1][2] * sys_in[2][1]
                  + mtx[2][2] * sys_in[2][2];
   
    // Copy sys_tmp to sys_out if rotating in place 
   if (sys_in == sys_out) {
      memcpy(sys_out, sys_tmp, sizeof(double)*9); 
   }
}
void AnalyticGeometryTool::rotate_system ( double  mtx[4][4],
double  sys_in[4][4],
double  sys_out[4][4] 
)

Definition at line 628 of file AnalyticGeometryTool.cpp.

{
   double sys_tmp[4][4];
   double* p_sys_tmp;
   
   // Check to see if rotating in place 
   if (sys_in == sys_out) {
      copy_mtx( sys_in, sys_tmp );
      p_sys_tmp = (double *)sys_tmp; 
   }
   else 
      // Have p_sys_tmp point at outgoing memory location 
      p_sys_tmp = (double *)sys_out;
   
   
   // X-vector 
   p_sys_tmp[0] =   mtx[0][0] * sys_in[0][0]
                  + mtx[1][0] * sys_in[0][1]
                  + mtx[2][0] * sys_in[0][2];
   p_sys_tmp[1] =   mtx[0][1] * sys_in[0][0]
                  + mtx[1][1] * sys_in[0][1]
                  + mtx[2][1] * sys_in[0][2];
   p_sys_tmp[2] =   mtx[0][2] * sys_in[0][0]
                  + mtx[1][2] * sys_in[0][1]
                  + mtx[2][2] * sys_in[0][2];   
   
   // Y-vector 
   p_sys_tmp[4] =   mtx[0][0] * sys_in[1][0]
                  + mtx[1][0] * sys_in[1][1]
                  + mtx[2][0] * sys_in[1][2];
   p_sys_tmp[5] =   mtx[0][1] * sys_in[1][0]
                  + mtx[1][1] * sys_in[1][1]
                  + mtx[2][1] * sys_in[1][2];
   p_sys_tmp[6] =   mtx[0][2] * sys_in[1][0]
                  + mtx[1][2] * sys_in[1][1]
                  + mtx[2][2] * sys_in[1][2];
   
   // Z-vector 
   p_sys_tmp[8] =   mtx[0][0] * sys_in[2][0]
                  + mtx[1][0] * sys_in[2][1]
                  + mtx[2][0] * sys_in[2][2];
   p_sys_tmp[9] =   mtx[0][1] * sys_in[2][0]
                  + mtx[1][1] * sys_in[2][1]
                  + mtx[2][1] * sys_in[2][2];
   p_sys_tmp[10] =  mtx[0][2] * sys_in[2][0]
                  + mtx[1][2] * sys_in[2][1]
                  + mtx[2][2] * sys_in[2][2];
   
   // Maintain the origin 
   p_sys_tmp[3] = sys_in[0][3];
   p_sys_tmp[7] = sys_in[1][3];
   p_sys_tmp[11] = sys_in[2][3];
   p_sys_tmp[15] = sys_in[3][3];   
   p_sys_tmp[12] = sys_in[3][0];
   p_sys_tmp[13] = sys_in[3][1];    
   p_sys_tmp[14] = sys_in[3][2];   
   
    // Copy sys_tmp to sys_out if rotating in place 
   if (sys_in == sys_out) {
      memcpy(sys_out, sys_tmp, sizeof(double)*16);
   }
}
void AnalyticGeometryTool::round_near_val ( double &  val)

Round value to either -1, 0, or 1 if within epsilon to to one of these. Epsilon determined by get_epsilon/set_epsilon functions.

Definition at line 157 of file AnalyticGeometryTool.cpp.

{
   if( dbl_eq( val, 0.0 ) )
      val = 0.0;
   else if( dbl_eq( val, 1.0 ) )
      val  = 1.0;
   else if( dbl_eq( val, -1.0 ) )
      val = -1.0;
}
double AnalyticGeometryTool::set_epsilon ( double  new_epsilon) [inline]

Set epsilon used for double number comparisons. Default value is 1e-6.

Definition at line 995 of file AnalyticGeometryTool.hpp.

{ 
   double old_epsilon = agtEpsilon;
   agtEpsilon = new_epsilon;
   return old_epsilon;
}
void AnalyticGeometryTool::set_pnt ( double  pnt[3],
double  x,
double  y,
double  z 
)

Sets the value of pnt to x,y,z (inline function).

void AnalyticGeometryTool::setup_arc ( double  vec_1[3],
double  vec_2[3],
double  origin[3],
double  start_angle,
double  end_angle,
double  radius,
AGT_3D_Arc arc 
)

Functions to populate an arc structure. The arc is defined with two orthogonal vectors in a plane, the arc origin, the beginning angle in radians to rotate from the start vector (towards second vector), the ending angle in radians to rotate from the start vector, and the radius of the arc. The arc is parameterized from 0 to 1.

Definition at line 1735 of file AnalyticGeometryTool.cpp.

{
   copy_pnt( vec_1, arc.Vec1 );
   copy_pnt( vec_2, arc.Vec2 );
   copy_pnt( origin, arc.Origin );
   arc.StartAngle = start_angle;
   arc.EndAngle = end_angle;
   arc.Radius = radius;
}
void AnalyticGeometryTool::setup_arc ( CubitVector vec_1,
CubitVector vec_2,
CubitVector  origin,
double  start_angle,
double  end_angle,
double  radius,
AGT_3D_Arc arc 
)

Functions to populate an arc structure. The arc is defined with two orthogonal vectors in a plane, the arc origin, the beginning angle in radians to rotate from the start vector (towards second vector), the ending angle in radians to rotate from the start vector, and the radius of the arc. The arc is parameterized from 0 to 1.

Definition at line 1748 of file AnalyticGeometryTool.cpp.

{
   vec_1.get_xyz( arc.Vec1 );
   vec_2.get_xyz( arc.Vec2 );
   origin.get_xyz( arc.Origin );
   arc.StartAngle = start_angle;
   arc.EndAngle = end_angle;
   arc.Radius = radius;
}
double AnalyticGeometryTool::Sign ( double  x) [private]
double AnalyticGeometryTool::Sin ( double  x) [private]

Definition at line 5988 of file AnalyticGeometryTool.cpp.

{
    double c[3];
    int positive = 0, ip[3];
    int negative = 0, in[3];
    int i;
    for (i = 0; i < 3; i++)
    {
        c[i] = line->N.x*T->v[i].x + line->N.y*T->v[i].y - line->c;
        if ( c[i] > 0.0 )
        {
            ip[positive++] = i;
        }
        else if ( c[i] < 0.0 )
        {
            in[negative++] = i;
        }
    }

    // For a split to occur, one of the c_i must be positive and one must
    // be negative.
    AgtTriList* list = NULL;

    if ( negative == 0 ) // T is completely on the positive side of line
        return 0;

    if ( positive == 0 )
    {
        // T is completely on the negative side of line
        list = new AgtTriList;
        list->tri = T;
        list->next = 0;
        return list;
    }

    // T is split by line.  Determine how it is split and how to decompose
    // the negative-side portion into triangles (3 cases).
    double w0, w1, cdiff;
    Point2 intrp[2];

    if ( positive == 2 )
    {
        // ++-
        for (i = 0; i < 2 /* = positive */; i++)
        {
            cdiff = c[ip[i]]-c[in[0]];
            w0 = -c[in[0]]/cdiff;
            w1 = c[ip[i]]/cdiff;
            T->v[ip[i]].x = w0*T->v[ip[i]].x+w1*T->v[in[0]].x;
            T->v[ip[i]].y = w0*T->v[ip[i]].y+w1*T->v[in[0]].y;
        }
        list = new AgtTriList;
        list->tri = T;
        list->next = 0;
    }
    else if ( positive == 1 )
    {
        if ( negative == 2 )
        {
            // +--
            for (i = 0; i < 2 /* = negative */; i++)
            {
                cdiff = c[ip[0]]-c[in[i]];
                w0 = -c[in[i]]/cdiff;
                w1 = c[ip[0]]/cdiff;
                intrp[i].x = w0*T->v[ip[0]].x+w1*T->v[in[i]].x;
                intrp[i].y = w0*T->v[ip[0]].y+w1*T->v[in[i]].y;
            }

            T->v[ip[0]] = intrp[0];
            list = new AgtTriList;
            list->tri = T;

            Triangle* T1 = new Triangle;
            T1->v[0] = intrp[0];
            T1->v[1] = T->v[in[1]];
            T1->v[2] = intrp[1];
            list->next = new AgtTriList;
            list->next->tri = T1;
            list->next->next = 0;
        }
        else
        {
            // +-0
            cdiff = c[ip[0]]-c[in[0]];
            w0 = -c[in[0]]/cdiff;
            w1 = c[ip[0]]/cdiff;
            T->v[ip[0]].x = w0*T->v[ip[0]].x+w1*T->v[in[0]].x;
            T->v[ip[0]].y = w0*T->v[ip[0]].y+w1*T->v[in[0]].y;
            list = new AgtTriList;
            list->tri = T;
            list->next = 0;
        }
    }

    return list;
}
double AnalyticGeometryTool::Sqrt ( double  x) [inline, private]

Definition at line 975 of file AnalyticGeometryTool.hpp.

{
    return double(sqrt(x));
}
void AnalyticGeometryTool::transform_line ( double  rot_mtx[4][4],
double  origin[3],
double  axis[3] 
)

Transforms a line using the given transformation mtx. The mtx is typically built using add_origin_to_rotation_mtx.

Definition at line 237 of file AnalyticGeometryTool.cpp.

{
   double end_pnt[3]; // Find arbitrary end point on line
   next_pnt( origin, axis, 10.0, end_pnt );
   
   transform_pnt( rot_mtx, origin, origin );
   transform_pnt( rot_mtx, end_pnt, end_pnt );
   
   axis[0] = end_pnt[0] - origin[0];
   axis[1] = end_pnt[1] - origin[1];
   axis[2] = end_pnt[2] - origin[2];
}
void AnalyticGeometryTool::transform_line ( double  rot_mtx[4][4],
CubitVector origin,
CubitVector axis 
)

Transforms a line using the given transformation mtx. The mtx is typically built using add_origin_to_rotation_mtx.

Definition at line 251 of file AnalyticGeometryTool.cpp.

{
   CubitVector end_point; // Find arbitrary end point on line
   origin.next_point( axis, 10.0, end_point );
   
   double start_pnt[3], end_pnt[3];
   copy_pnt( origin, start_pnt );
   copy_pnt( end_point, end_pnt );
   
   transform_pnt( rot_mtx, start_pnt, start_pnt );
   transform_pnt( rot_mtx, end_pnt, end_pnt );
   
   axis.x( end_pnt[0] - start_pnt[0] );
   axis.y( end_pnt[1] - start_pnt[1] );
   axis.z( end_pnt[2] - start_pnt[2] );
   
   origin.set( start_pnt );
}
void AnalyticGeometryTool::transform_pnt ( double  m[4][4],
double  pin[3],
double  pout[3] 
)

This functions applies the transformation matrix to the specified point and returns the new point coordinates through the call list. Point in and point out can be the same memory address or different.

Definition at line 170 of file AnalyticGeometryTool.cpp.

{
    double p[3];  // working buffer
    
    // Check if transformation can occur 
    if (!m) {
       if (pin && pout)
      copy_pnt(pin, pout);
       return;
    }
    
    // Perform transformation    
    p[0] = m[0][0] * pin[0] + m[1][0] * pin[1] + m[2][0] * pin[2] + m[3][0];
    p[1] = m[0][1] * pin[0] + m[1][1] * pin[1] + m[2][1] * pin[2] + m[3][1];
    p[2] = m[0][2] * pin[0] + m[1][2] * pin[1] + m[2][2] * pin[2] + m[3][2];
    
    // Copy work buffer to out point    
    copy_pnt(p,pout);
} 
void AnalyticGeometryTool::transform_vec ( double  m3[3][3],
double  vin[3],
double  vout[3] 
)

This routine applies a transformation matrix to a specified vector and returns the new vector orientation. Vector in and vector out can be the same memory address or different.

Definition at line 192 of file AnalyticGeometryTool.cpp.

{
   double v[3];    // working buffer
   
   // Determine if transformation can occur 
   if (!m3) {      
      if (vin && vout)
         copy_vec(vin, vout);
      return;
   }
   
   // Perform transformation 
   v[0] = m3[0][0] * vin[0] + m3[1][0] * vin[1] + m3[2][0] * vin[2];
   v[1] = m3[0][1] * vin[0] + m3[1][1] * vin[1] + m3[2][1] * vin[2];
   v[2] = m3[0][2] * vin[0] + m3[1][2] * vin[1] + m3[2][2] * vin[2];
   
   // Copy work buffer to vector out    
   copy_pnt(v,vout);
}
void AnalyticGeometryTool::transform_vec ( double  m4[4][4],
double  vin[3],
double  vout[3] 
)

This routine applies a transformation matrix to a specified vector and returns the new vector orientation. Vector in and vector out can be the same memory address or different.

Definition at line 214 of file AnalyticGeometryTool.cpp.

{
   double v[3];    // working buffer
   
   // Determine if transformation can occur 
   if (!m4) {      
      if (vin && vout)
         copy_vec(vin, vout);
      return;
   }
   
   // Perform transformation
   
   v[0] = m4[0][0] * vin[0] + m4[1][0] * vin[1] + m4[2][0] * vin[2];
   v[1] = m4[0][1] * vin[0] + m4[1][1] * vin[1] + m4[2][1] * vin[2];
   v[2] = m4[0][2] * vin[0] + m4[1][2] * vin[1] + m4[2][2] * vin[2];
   
   // Copy work buffer to vector out    
   copy_pnt(v,vout);
}
void AnalyticGeometryTool::TriangleLines ( Triangle T,
AgtLine  line[3] 
) [private]

Definition at line 5957 of file AnalyticGeometryTool.cpp.

{
    line[0] = EdgeToLine(&T->v[0],&T->v[1]);
    line[1] = EdgeToLine(&T->v[0],&T->v[2]);
    line[2] = EdgeToLine(&T->v[1],&T->v[2]);

    // make sure normals point outwards
    Point2 avr = { 0.0, 0.0 };
    int i;
    for (i = 0; i < 3; i++)
    {
        avr.x += T->v[i].x;
        avr.y += T->v[i].y;
    }
    static double oneThird = 1.0/3.0;
    avr.x *= oneThird;
    avr.y *= oneThird;

    for (i = 0; i < 3; i++)
    {
        double dot = line[i].N.x*avr.x+line[i].N.y*avr.y-line[i].c;
        if ( dot > 0.0 )
        {
            line[i].N.x = -line[i].N.x;
            line[i].N.y = -line[i].N.y;
            line[i].c = -line[i].c;
        }
    }
}
CubitStatus AnalyticGeometryTool::unit_vec ( double  vin[3],
double  vout[3] 
)

Finds unit vector of input vector (vout can equal vin - converts in place).

Definition at line 1093 of file AnalyticGeometryTool.cpp.

{
   double rmag; // holds magnitude of vector 
      
   // Calculate vector magnitude 
   rmag = sqrt(vin[0]*vin[0] + vin[1]*vin[1] + vin[2]*vin[2]);
   
   // check if vector has a magnitude - catch division by zero 
   
   if (dbl_eq(rmag, 0.0)) {
      if (vin != vout)
     copy_pnt(vin, vout);
      return CUBIT_FAILURE;
   }
   
   // divide each element of the vector by the magnitude 
      
   vout[0] = vin[0] / rmag;
   vout[1] = vin[1] / rmag;
   vout[2] = vin[2] / rmag;
      
   return CUBIT_SUCCESS;
}
double AnalyticGeometryTool::UnitRandom ( ) [inline, private]

Definition at line 980 of file AnalyticGeometryTool.hpp.

{
    return double(rand())/double(RAND_MAX);
}
void AnalyticGeometryTool::vecs_to_mtx ( double  xvec[3],
double  yvec[3],
double  zvec[3],
double  matrix[3][3] 
)

Creates a transformation matrix given three vectors (and origin for 4x4 case).

Definition at line 923 of file AnalyticGeometryTool.cpp.

{      
   if (xvec) 
      copy_pnt(xvec, matrix[0]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[0], matrix[0]);
   
   if (yvec)
      copy_pnt(yvec, matrix[1]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[1], matrix[1]);
   
   if (zvec)
      copy_pnt(zvec, matrix[2]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[2], matrix[2]);
}   
void AnalyticGeometryTool::vecs_to_mtx ( double  xvec[3],
double  yvec[3],
double  zvec[3],
double  origin[3],
double  matrix[4][4] 
)

Creates a transformation matrix given three vectors (and origin for 4x4 case).

Definition at line 944 of file AnalyticGeometryTool.cpp.

{      
   if (xvec) 
      copy_pnt(xvec, matrix[0]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[0], matrix[0]);
   
   if (yvec)
      copy_pnt(yvec, matrix[1]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[1], matrix[1]);
   
   if (zvec)
      copy_pnt(zvec, matrix[2]);
   else
      copy_pnt(AGT_IDENTITY_MTX_3X3[2], matrix[2]);
   
   if( origin )
      copy_pnt(origin, matrix[3]);
   else
   {
      matrix[3][0] = 0.0;
      matrix[3][1] = 0.0;
      matrix[3][2] = 0.0;
   }

   matrix[0][3] = 0.0;
   matrix[1][3] = 0.0;
   matrix[2][3] = 0.0;
   matrix[3][3] = 1.0;    
}
double AnalyticGeometryTool::Volume ( int  N,
Point3 pt,
double  angle[3] 
) [private]

Definition at line 2845 of file AnalyticGeometryTool.cpp.

{
    double cs0 = cos(angle[0]), sn0 = sin(angle[0]);
    double cs1 = cos(angle[1]), sn1 = sin(angle[1]);
    double axis[3] = { cs0*sn1, sn0*sn1, cs1 };
    double rot[3][3];
    AngleAxisToMatrix(angle[2],axis,rot);

    double min[3] =
    {
        rot[0][0]*pt[0].x+rot[1][0]*pt[0].y+rot[2][0]*pt[0].z,
        rot[0][1]*pt[0].x+rot[1][1]*pt[0].y+rot[2][1]*pt[0].z,
        rot[0][2]*pt[0].x+rot[1][2]*pt[0].y+rot[2][2]*pt[0].z
    };

    double max[3] = { min[0], min[1], min[2] };

    for (int i = 1; i < N; i++)
    {
        double test[3] =
        {
            rot[0][0]*pt[i].x+rot[1][0]*pt[i].y+rot[2][0]*pt[i].z,
            rot[0][1]*pt[i].x+rot[1][1]*pt[i].y+rot[2][1]*pt[i].z,
            rot[0][2]*pt[i].x+rot[1][2]*pt[i].y+rot[2][2]*pt[i].z
        };

        if ( test[0] < min[0] )
            min[0] = test[0];
        else if ( test[0] > max[0] )
            max[0] = test[0];

        if ( test[1] < min[1] )
            min[1] = test[1];
        else if ( test[1] > max[1] )
            max[1] = test[1];

        if ( test[2] < min[2] )
            min[2] = test[2];
        else if ( test[2] > max[2] )
            max[2] = test[2];
    }

    double volume = (max[0]-min[0])*(max[1]-min[1])*(max[2]-min[2]);
    return volume;
}

Member Data Documentation

Definition at line 760 of file AnalyticGeometryTool.hpp.

Definition at line 758 of file AnalyticGeometryTool.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines