C++ CubitVector::z方法代码示例

 //=========================================================================== 
//Function Name: intersect 
// 
//Member Type:  PUBLIC 
//Description:  intersect the edge with a segment.  Assumes segment and edge
//              are on the same plane (project to facet plane first) 
//=========================================================================== 
CubitStatus CubitFacetEdge::intersect(
  CubitVector &aa, CubitVector &bb, // end point of segment
  CubitVector &norm,  // normal of the common plane
  CubitVector &qq,  // return the intersection point 
  CubitBoolean &does_intersect ) // return status of intersection
{  
 
  CubitPoint *pt0 = this->point(0); 
  CubitPoint *pt1 = this->point(1); 

  double P0[2], P1[2], AA[2], BB[2];
  CubitVector absnorm(fabs(norm.x()), fabs(norm.y()), fabs(norm.z()));
  if (absnorm.x() >= absnorm.y() && absnorm.x() >= absnorm.z())
  {
    P0[0] = pt0->coordinates().y();  P0[1] = pt0->coordinates().z();
    P1[0] = pt1->coordinates().y();  P1[1] = pt1->coordinates().z();
    AA[0] = aa.y();                  AA[1] = aa.z();
    BB[0] = bb.y();                  BB[1] = bb.z();
  }
  else if (absnorm.y() >= absnorm.x() && absnorm.y() >= absnorm.z())
  {
    P0[0] = pt0->coordinates().z();  P0[1] = pt0->coordinates().x();
    P1[0] = pt1->coordinates().z();  P1[1] = pt1->coordinates().x();
    AA[0] = aa.z();                  AA[1] = aa.x();
    BB[0] = bb.z();                  BB[1] = bb.x();
  }
  else
  {
    P0[0] = pt0->coordinates().x();  P0[1] = pt0->coordinates().y();
    P1[0] = pt1->coordinates().x();  P1[1] = pt1->coordinates().y();
    AA[0] = aa.x();                  AA[1] = aa.y();
    BB[0] = bb.x();                  BB[1] = bb.y();
  }

  double QQ[4], s;
  int ninter = intersect_2D_segments(P0, P1, AA, BB, QQ);

  if (ninter != 1)
  {
    does_intersect = CUBIT_FALSE;
    return CUBIT_SUCCESS;
  }
  does_intersect = CUBIT_TRUE;
 
  double dx = P1[0] - P0[0];
  double dy = P1[1] - P0[1];
  if (fabs(dx) > fabs(dy))
    s = (QQ[0] - P0[0]) / dx;
  else
    s = (QQ[1] - P0[1]) / dy;
  
  qq = pt0->coordinates() + s * (pt1->coordinates() - pt0->coordinates());
 
  return CUBIT_SUCCESS; 
}  

template <class Z> KDDTreeNode<Z> *KDDTree<Z>::recursive_balance
 (DIMENSION dim, int left, int right, KDDTreeNode<Z>** array, KDDTreeNode<Z>* parent)
{
  if (left > right)
  {
    return NULL;
  }
  else
  {
    KDDTreeNode<Z>* P;
    int K;

    if (left != right)
    {
      K = modifind (dim, left, right, array);
      P = array[K];
    }
    else
    {
      K = left;
      P = array[left];
    }

    myNodeList.push (P);

    P->safetyBox.reset (P->boundingBox);
    for (int i = left; i <= right; i++)
    {
      CubitVector imin = array[i]->safetyBox.minimum();
      CubitVector imax = array[i]->safetyBox.maximum();
      CubitVector pmin = P->safetyBox.minimum();
      CubitVector pmax = P->safetyBox.maximum();

      if (imin.x() < pmin.x()) pmin.x (imin.x());
      if (imin.y() < pmin.y()) pmin.y (imin.y());
      if (imin.z() < pmin.z()) pmin.z (imin.z());
      if (imax.x() > pmax.x()) pmax.x (imax.x());
      if (imax.y() > pmax.y()) pmax.y (imax.y());
      if (imax.z() > pmax.z()) pmax.z (imax.z());

      P->safetyBox.reset (pmin, pmax);
    }

    DIMENSION nextDim;
    switch (dim)
    {
      case DIMX: nextDim = DIMY; break;
      case DIMY: nextDim = DIMZ; break;
      default:   nextDim = DIMX;
    }

    P->set_disc (dim);
    P->parent = parent;

    P->left = recursive_balance (nextDim, left, K - 1, array, P);
    P->right = recursive_balance (nextDim, K + 1, right, array, P);

    return P;
  }
}

//-------------------------------------------------------------------------
// Purpose       : Compute and return the bounding box for this sphere.
//
// Special Notes :
//
//-------------------------------------------------------------------------
CubitBox SphereEvaluator::bounding_box( void ) const
{
    CubitVector center = mTmatrix * mEvalData.center;
            
    CubitVector min( center.x() - mEvalData.radius,
                     center.y() - mEvalData.radius,
                     center.z() - mEvalData.radius );
    CubitVector max( center.x() + mEvalData.radius,
                     center.y() + mEvalData.radius,
                     center.z() + mEvalData.radius );

    CubitBox thebbox( min, max );

    return thebbox;
}

CubitPointContainment OCCSurface::point_containment( const CubitVector &point )
{
   TopoDS_Face *face = get_TopoDS_Face();
   gp_Pnt p(point.x(), point.y(), point.z());
   double tol = OCCQueryEngine::instance()->get_sme_resabs_tolerance();

   //It's checking the state of the projected point of THIS Point
   BRepClass_FaceClassifier face_classifier;
   face_classifier.Perform(*face, p, tol);
   TopAbs_State state = face_classifier.State();
   
   //if surface is part of a periodic TopoDS_Face, it'll check the point
   //againt the whole periodic Face, even it outside the occsurface 
   //boundary, if it's on its periodic extension, it'll return as in. 
   if (state == TopAbs_IN)
   {
     //double check if the point is projected on the surface
     CubitVector closest_point;
     this->closest_point_trimmed(point, closest_point);
     if(point.distance_between(closest_point) < tol) 
       return CUBIT_PNT_INSIDE;
     else
       return CUBIT_PNT_OUTSIDE;
   }
   else if (state == TopAbs_OUT)
     return CUBIT_PNT_OUTSIDE;
   else if (state == TopAbs_ON)
     return CUBIT_PNT_BOUNDARY;

   return CUBIT_PNT_UNKNOWN;
}

CubitStatus OCCSurface::principal_curvatures(
  CubitVector const& location, 
  double& curvature_1,
  double& curvature_2,
  CubitVector* closest_location )
{
  BRepAdaptor_Surface asurface(*myTopoDSFace);
  gp_Pnt p(location.x(), location.y(), location.z()), newP(0.0, 0.0, 0.0);
  double minDist=0.0, u, v;
  int i;
  BRepLProp_SLProps SLP(asurface, 2, Precision::PConfusion());
  Extrema_ExtPS ext(p, asurface, Precision::Approximation(), Precision::Approximation());
  if (ext.IsDone() && (ext.NbExt() > 0)) {
	  for ( i = 1 ; i <= ext.NbExt() ; i++ ) {
		  if ( (i==1) || (p.Distance(ext.Point(i).Value()) < minDist) ) {
			  minDist = p.Distance(ext.Point(i).Value());
			  newP = ext.Point(i).Value();
			  ext.Point(i).Parameter(u, v);
			  SLP.SetParameters(u, v);
		  }
	  }
  }
  if (closest_location != NULL)
    *closest_location = CubitVector(newP.X(), newP.Y(), newP.Z());

  if (SLP.IsCurvatureDefined())
  {
    curvature_1 = SLP.MinCurvature();
    curvature_2 = SLP.MaxCurvature();
  }
  return CUBIT_SUCCESS;
}

double SurfaceOverlapFacet::distance_from_position( CubitVector &position )
{
  double s,t;  
  Point3 tmp_point;
  tmp_point.x = position.x();
  tmp_point.y = position.y();
  tmp_point.z = position.z();
  return agt->MinPointTriangle( tmp_point, this->t, s, t );
} 

//--------------------------------------------------------------------------
//Algorithm: min_dist_sq
//Description:  Finds the minimum distance squared between the given
//              point and the box. If the point is on or in the box, the
//              min distance is zero.
//--------------------------------------------------------------------------
template <class Z> MY_INLINE
double RStarTree<Z>::min_dist_sq(CubitVector &q,
                             CubitBox &b_box)
{
  CubitVector b_min, b_max;
  b_min = b_box.minimum();
  b_max = b_box.maximum();
  double dist;
  CubitVector r;

  if ( q.x() < b_min.x() )
    r.x(b_min.x());
  else if ( q.x() > b_max.x() )
    r.x(b_max.x());
  else
    r.x(q.x());
  
  if ( q.y() < b_min.y() )
    r.y(b_min.y());
  else if ( q.y() > b_max.y() )
    r.y(b_max.y());
  else
    r.y(q.y());
  
  if ( q.z() < b_min.z() )
    r.z(b_min.z());
  else if ( q.z() > b_max.z() )
    r.z(b_max.z());
  else
    r.z(q.z());
  
  dist = (q-r).length_squared();

  return dist;
}

//=========================================================================== 
//Function Name: closest_point 
// 
//Member Type:  PUBLIC 
//Description:  return the closest point on segment defined by the edge 
//=========================================================================== 
CubitStatus CubitFacetEdge::closest_point(const CubitVector &point,  
                                          CubitVector &closest_point ) 
{  
  //CubitStatus rv = CUBIT_SUCCESS; 
  CubitPoint *pt0 = this->point(0); 
  CubitPoint *pt1 = this->point(1); 
 
  // the edge vector 
 
  CubitVector e0 ( pt1->x() - pt0->x(), 
                   pt1->y() - pt0->y(), 
                   pt1->z() - pt0->z() ); 
  double elen = e0.normalize(); 
   
  // vector from vert0 to point 
 
  CubitVector v0 ( point.x() - pt0->x(), 
                   point.y() - pt0->y(), 
                   point.z() - pt0->z() ); 
   
  // project to edge 
 
  double len = v0%e0; 
  if (len <= 0.0) 
  { 
    closest_point = pt0->coordinates(); 
  } 
  else if( len >= elen ) 
  { 
    closest_point = pt1->coordinates(); 
  } 
  else 
  { 
    closest_point.x ( pt0->x() + e0.x() * len ); 
    closest_point.y ( pt0->y() + e0.y() * len ); 
    closest_point.z ( pt0->z() + e0.z() * len ); 
  } 
 
  return CUBIT_SUCCESS; 
}

//-------------------------------------------------------------------------
// Purpose       : make arbitrary parameterization
//
// Special Notes : 
//
// Creator       : Jason Kraftcheck 
//
// Creation Date : 07/06/98
//-------------------------------------------------------------------------
void ParamCubitPlane::make_parameterization()
{
	//Choose the zero-point for the parameterization
	//as close to the origin as possible.
//	p_.set( 0.0, 0.0, 0.0);
//	move_to_plane( p_ );
	is_plane_valid_ = CUBIT_TRUE;
	
	const CubitVector temp_p(0.0, 0.0, 0.0);
	CubitStatus s = closest_point( temp_p, p_ );
	assert( s == CUBIT_SUCCESS );
	
	CubitVector n = normal();
	CubitVector p1;
	
	p1 = p_;
	double x = fabs( n.x() );
	double y = fabs( n.y() );
	double z = fabs( n.z() );
	
	//Choose a direction from the zero point (p_) for
	//the second point as the direction of the smallest
	//component of the normal.  The third point defining
	//the plane will be defined by the cross product of
	//the vector from the zero_point to this point and
	//the normal vector of the plane.
	if( (x <= y) && (x <= z) )
	{
		p1.x( p1.x() + 1 );
	}
	else if( (y <= x) && (y <= z) )
	{
		p1.y( p1.y() + 1 );
	}
	else
	{
		p1.z( p1.z() + 1 );
	}
	
	move_to_plane( p1 );
	s_ = p1 - p_;
	t_ = -(s_ * n);
	n_ = s_ * t_;
	
	double n_len = n_.length();
	if( 1000 * CUBIT_DBL_MIN > n_len ) n_epsilon_ = CUBIT_DBL_MIN;
	else n_epsilon_ = n_len / 1000;
	
 
 	is_plane_valid_= CUBIT_TRUE;
}

//===================================================================================
// Function: transform_to_uv_local (Private)
// Description: It transforms points from the Best Fit plane plane to UV space
// Author: Shiraj Khan
// Date: 1/21/2003
//===================================================================================
CubitStatus LoopParamTool::transform_to_uv_local(CubitVector &xyz_location, CubitVector &uv_location) 
{
    // Translate to local origin at center

  CubitVector vect = xyz_location - uvCenter;

    // Multiply by transpose (inverse) of transformation vector

  uv_location.x( vect % Du );
  uv_location.y( vect % Dv );
  uv_location.z( 1.0 );

  return CUBIT_SUCCESS;
}

//-------------------------------------------------------------------------
// Purpose       : Computes the closest_point on the surface to the input 
//                 location.  Optionally, it also computes and returns
//                 the normal to the surface at closest_location and the 
//                 principal curvatures(1-min, 2-max)
//
//-------------------------------------------------------------------------
CubitStatus OCCSurface::closest_point( CubitVector const& location, 
                                         CubitVector* closest_location,
                                         CubitVector* unit_normal_ptr,
                                         CubitVector* curvature_1,
                                         CubitVector* curvature_2)
{
  BRepAdaptor_Surface asurface(*myTopoDSFace);
  gp_Pnt p(location.x(), location.y(), location.z()), newP(0.0, 0.0, 0.0);
  double minDist=0.0, u, v;
  int i;
  BRepLProp_SLProps SLP(asurface, 2, Precision::PConfusion());
  Extrema_ExtPS ext(p, asurface, Precision::Approximation(), Precision::Approximation());
  if (ext.IsDone() && (ext.NbExt() > 0)) {
	  for ( i = 1 ; i <= ext.NbExt() ; i++ ) {
		  if ( (i==1) || (p.Distance(ext.Point(i).Value()) < minDist) ) {
			  minDist = p.Distance(ext.Point(i).Value());
			  newP = ext.Point(i).Value();
			  ext.Point(i).Parameter(u, v);
			  SLP.SetParameters(u, v);
		  }
	  }
  
	if (closest_location != NULL)
 	 	*closest_location = CubitVector(newP.X(), newP.Y(), newP.Z());
  	if (unit_normal_ptr != NULL) {
	  gp_Dir normal;
          //normal of a RefFace point to outside of the material
	  if (SLP.IsNormalDefined()) {
	    normal = SLP.Normal();
            CubitSense sense = get_geometry_sense();
            if(sense == CUBIT_REVERSED)
              normal.Reverse() ;
	      *unit_normal_ptr = CubitVector(normal.X(), normal.Y(), normal.Z()); 
	  }
  	}
  
        gp_Dir MaxD, MinD;
        if (SLP.IsCurvatureDefined())
        {
	   SLP.CurvatureDirections(MaxD, MinD);
           if (curvature_1 != NULL)
              *curvature_1 = CubitVector(MinD.X(), MinD.Y(), MinD.Z());
           if (curvature_2 != NULL)
              *curvature_2 = CubitVector(MaxD.X(), MaxD.Y(), MaxD.Z());
        }
  }
  return CUBIT_SUCCESS;
}

//- "insert_node"
//- Dynamically insert the data into the k-d tree 
//- 
//- Algorithm INSERT (From Bentley):
//-   This algoritm is passed an object "data" of class "Z",
//-   which has a bounding_box() method.  If there is already
//-   a node in the tree with equal bounding box center point,
//-   it is put in the right subtree.
//-   I0. [Create new node] Create a node P with the bounding box
//-       specified, and set P.LEFT <- null, P.RIGHT <- null, and
//-       P.DISC <- null.
//-   I1. [Check for null tree] If ROOT = null, then set ROOT <- P
//-       and return CUBIT_SUCCESS; otherwise, set Q <- ROOT (Q
//-       will move down the tree).
//-   I2. [Compare] Compare the nodes and set the child in the
//-       correct direction to T.
//-   I3. [Move down] Set Q <- child of Q and go to I2.
//-   I4. [Insert new node in tree] Set the child of Q to P, then
//-       set the children of P to null. Set the discriminator of
//-       P to be the discriminator after that in Q.
//-
template <class Z> CubitStatus KDDTree<Z>::insert_node (KDDTreeNode<Z>* P)
{
  KDDTreeNode<Z> *F = NULL;      // father node
  KDDTreeNode<Z> *T;             // temp node

  if (root == NULL)
  {
    root = P;
    P->set_disc (DIMX);
  }
  else
  {
    T = root;
    DIRECTION direction = DIR_NULL;

    while (T != NULL)
    {
      F = T; // remember the father
      direction = P->compare (T);

      CubitVector tmin = T->safetyBox.minimum();
      CubitVector tmax = T->safetyBox.maximum();
      CubitVector pmin = P->safetyBox.minimum();
      CubitVector pmax = P->safetyBox.maximum();

      if (pmin.x() < tmin.x()) tmin.x (pmin.x());
      if (pmin.y() < tmin.y()) tmin.y (pmin.y());
      if (pmin.z() < tmin.z()) tmin.z (pmin.z());
      if (pmax.x() > tmax.x()) tmax.x (pmax.x());
      if (pmax.y() > tmax.y()) tmax.y (pmax.y());
      if (pmax.z() > tmax.z()) tmax.z (pmax.z());

      T->safetyBox.reset (tmin, tmax);

      T = T->get_child (direction);
    }

    F->set_child (P, direction);
    P->set_disc (F->next_disc());
    P->parent = F;
  }
  myNodeList.push (P);

  return CUBIT_SUCCESS;      
}

//===================================================================================
// Function: transform_to_xyz_local (Priavate)
// Description: same as title
// Author: Shiraj Khan
// Date: 1/21/2003
//===================================================================================
CubitStatus LoopParamTool::transform_to_xyz_local(CubitVector &xyz_location, CubitVector &uv_location) 
{
// Multiply by transformation matrix

  CubitVector vect;
  vect.x( uv_location.x() * Du.x() +
          uv_location.y() * Dv.x() );
  vect.y( uv_location.x() * Du.y() +
          uv_location.y() * Dv.y() );
  vect.z( uv_location.x() * Du.z() +
          uv_location.y() * Dv.z() );

    // Translate from origin

  xyz_location = vect + uvCenter;

  return CUBIT_SUCCESS;
}

//- Finds the minimum distance squared between the given point and the box. If
//-  the point is on or in the box, the min distance is zero.
template <class Z> double KDDTree<Z>::min_dist_sq (CubitVector &q, CubitBox &b_box)
{
  CubitVector b_min = b_box.minimum();
  CubitVector b_max = b_box.maximum();
  CubitVector r;

  //// set "r" in the x-dim
  if (q.x () < b_min.x ())
  {
    r.x (b_min.x ());
  }
  else if (q.x () > b_max.x ())
  {
    r.x (b_max.x ());
  }
  else
  {
    r.x (q.x ());
  }
  
  //// set "r" in the y-dim
  if (q.y () < b_min.y ())
  {
    r.y (b_min.y ());
  }
  else if (q.y () > b_max.y ())
  {
    r.y (b_max.y ());
  }
  else
  {
    r.y (q.y ());
  }
  
  //// set "r" in the z-dim
  if (q.z () < b_min.z ())
  {
    r.z (b_min.z ());
  }
  else if (q.z () > b_max.z ())
  {
    r.z (b_max.z ());
  }
  else
  {
    r.z (q.z ());
  }
  
  double dist = (q-r).length_squared();

  return dist;
}

//-------------------------------------------------------------------------
// Purpose       : Computes the closest_point on the trimmed surface to the 
//                 input location. 
//
// Special Notes : 
//-------------------------------------------------------------------------
void OCCSurface::closest_point_trimmed( CubitVector from_point, 
                                         CubitVector& point_on_surface)
{
  BRepAdaptor_Surface asurface(*myTopoDSFace);
  gp_Pnt p(from_point.x(), from_point.y(), from_point.z()), newP(0.0, 0.0, 0.0);
  double minDist=0.0;
  int i;
  Extrema_ExtPS ext(p, asurface, Precision::Approximation(), Precision::Approximation());
  if (ext.IsDone() && (ext.NbExt() > 0)) {
	  for ( i = 1 ; i <= ext.NbExt() ; i++ ) {
		  if ( (i==1) || (p.Distance(ext.Point(i).Value()) < minDist) ) {
			  minDist = p.Distance(ext.Point(i).Value());
			  newP = ext.Point(i).Value();
		  }
	  }
  }
  point_on_surface = CubitVector(newP.X(), newP.Y(), newP.Z());
}