C++ CubitVector类代码示例

//===========================================================================
//  Function: evaluate
//  Purpose:  evaluate the quad at a specified u,v coordinate
//  Notes:    Uses the interpolation order previously set up for the
//            facet-based surface.  The u,v system is based on the point
//            order and is defined as follows:
//
//            [3] =  (-1,1)  *============*  [2] = (1,1)
//                           |            |  
//                           |            | 
//                           |            |
//                           |            |
//                           |            |
//            [0] = (-1,-1)  *============*  [1] = (1,-1)
//
//  Date:     4/11/01
//  Author:   sjowen
//===========================================================================
CubitStatus CubitQuadFacet::evaluate( double u, double v, 
                                      CubitVector *eval_point,
                                      CubitVector *eval_normal,
                                      CubitVector *eval_du,
                                      CubitVector *eval_dv)
{
  CubitVector normal;
  CubitStatus stat = CUBIT_SUCCESS;
  int tri_index;
  double A, B, C;
  map_to_tri_system( u, v, tri_index, A, B, C );
  CubitVector ac(A,B,C);
  CubitFacet *tri_facet = get_tri_facet( tri_index );
  CubitVector point;
  stat = tri_facet->evaluate( ac, &point, eval_normal );
  if (!stat)
    return stat;
  if (eval_point != NULL)
    *eval_point = point;
  if (eval_du != NULL || eval_dv != NULL)
  {
    CubitVector du, dv;
    stat = eval_derivatives( u, v, point, du, dv );
    normal = du * dv;
    normal.normalize();
    if (eval_du != NULL)
      *eval_du = du;
    if (eval_dv != NULL)
      *eval_dv = dv;
    if (eval_normal != NULL)
      *eval_normal = normal;
  }
  return stat;
}

// Constructor
SurfaceOverlapFacet::SurfaceOverlapFacet( GPoint p[3] )
{
  t.b.x = p[0].x;
  t.b.y = p[0].y;
  t.b.z = p[0].z;

  t.e0.x = p[1].x - p[0].x;
  t.e0.y = p[1].y - p[0].y;
  t.e0.z = p[1].z - p[0].z;

  t.e1.x = p[2].x - p[0].x;
  t.e1.y = p[2].y - p[0].y;
  t.e1.z = p[2].z - p[0].z;

  CubitVector min;
   
  min.set( CUBIT_MIN_3( p[0].x, p[1].x, p[2].x ),
    CUBIT_MIN_3( p[0].y, p[1].y, p[2].y ),
    CUBIT_MIN_3( p[0].z, p[1].z, p[2].z ) );
  
  CubitVector max;
  
  max.set( CUBIT_MAX_3( p[0].x, p[1].x, p[2].x ),
    CUBIT_MAX_3( p[0].y, p[1].y, p[2].y ),
    CUBIT_MAX_3( p[0].z, p[1].z, p[2].z ) );

  boundingBox.reset( min, max );
}

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;
}

//-------------------------------------------------------------------------
// Purpose       : Computes the closest_point on the surface to the input 
//                 location.  Optionally, it also computes and returns
//                 the normal to the surface and the principal curvatures
//                 at closest_location.
//
//-------------------------------------------------------------------------
CubitStatus SphereEvaluator::closest_point( CubitVector const& location, 
                                            CubitVector* closest_location,
                                            CubitVector* unit_normal_ptr,
                                            CubitVector* curvature1_ptr,
                                            CubitVector* curvature2_ptr) const
{
    CubitTransformMatrix inverse_Tmatrix = mTmatrix;
    inverse_Tmatrix.inverse();
    CubitVector new_pt = inverse_Tmatrix * location;

    // ******************************************************
    // If requested, compute the normal at the input point.
    // ******************************************************

    if ( unit_normal_ptr != NULL )
    {
        CubitVector origin( 0.0, 0.0, 0.0 );
        CubitVector endpt = new_pt - mEvalData.center;

        origin = mTmatrix * origin;
        endpt  = mTmatrix * endpt;

        *unit_normal_ptr = endpt - origin;
        unit_normal_ptr->normalize();
    }

    // *************************************************************
    // If requested, compute the closest point to the input point.
    // *************************************************************

    if ( closest_location != NULL )
    {
        if ( location.within_tolerance( mEvalData.center, GEOMETRY_RESABS ) )
        {
            closest_location->set( mEvalData.center.x() + mEvalData.radius,
                                   mEvalData.center.y(),
                                   mEvalData.center.z() );
        }
        else
        {
            CubitVector vec = new_pt - mEvalData.center;
            vec.normalize();
            *closest_location = mEvalData.center + ( vec * mEvalData.radius );
            *closest_location = mTmatrix * (*closest_location);
        }
    }

    // ***********************************************************************
    // If requested, compute the curvature directions.
    // ***********************************************************************

    if ( curvature1_ptr &&
         curvature2_ptr )
    {
        PRINT_ERROR("Sphere's do not current support curvature pointer requests.\n");
        return CUBIT_FAILURE;
    }
  
    return CUBIT_SUCCESS;
}

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;
}

double
SurfaceOverlapFacet::projected_overlap( SurfaceOverlapFacet &other_facet, CubitBoolean draw_overlap )
{
  double tmp_double = agt->ProjectedOverlap( t, other_facet.t, draw_overlap );

  if( tmp_double > 0.00 ) 
  {
    CubitVector edge0(t.e0.x, t.e0.y, t.e0.z);  
    CubitVector edge1(t.e1.x, t.e1.y, t.e1.z);  
    CubitVector normal = edge0 * edge1;
    double area_facet1 = normal.length() / 2;

    edge0.set(other_facet.t.e0.x, other_facet.t.e0.y, other_facet.t.e0.z);  
    edge1.set(other_facet.t.e1.x, other_facet.t.e1.y, other_facet.t.e1.z);  
    normal = edge0 * edge1;
    double area_facet2 = normal.length() / 2;
    
    //don't report overlapping area between facets unless it is greater 
    //than one hundredth of the area of the smaller facet
    if( area_facet1 < area_facet2 )
    {
      if( tmp_double < (area_facet1*0.01))
        tmp_double = 0.0;
    }
    else if( tmp_double < (area_facet2*0.01 ))
      tmp_double = 0.0;
  }
  return tmp_double;
}

//
// Constructor / Destructor ...................................................
//
PointGridSearch::PointGridSearch( DLIList<CubitPoint*> *point_list,
                                  DLIList<CubitFacet*> *facet_list,
                                  float grid_scale)
{
    // find the geometric bounding range mesh
  assert( point_list && point_list->size() );
  bounding_range( *point_list );

    //find an estimate of the cell size based on the average square distance
    //of the edges on the surface of the geometry

  double grid_cell_size = 1.;
  
    // edge lengths
  double sum = 0.0;
  int i;
  maxEdgeLength = CUBIT_DBL_MIN;
  assert( facet_list || facet_list->size());
  for ( i = 0; i < facet_list->size(); i++)
  {
    CubitPoint *p1, *p2;
    facet_list->next(i)->get_edge_1(p1, p2);
    CubitVector temp = p2->coordinates() - p1->coordinates();
    double tmp_lsqrd =  temp.length_squared();
    sum += tmp_lsqrd;
  }
  grid_cell_size = sqrt( sum / facet_list->size() );

  maxEdgeLength = grid_cell_size;
    // evaluate grid parameters
  CubitVector cell(1.0, 1.0, 1.0);
  double grid_cell_width = grid_scale * grid_cell_size;
  cell *= ( GRID_EXTENSION + 5.0 ) * grid_cell_width;
  gridRangeMinimum = boundingRangeMinimum - cell;
  gridRangeMaximum = boundingRangeMaximum + cell;
  
  CubitVector range_width = gridRangeMaximum - gridRangeMinimum;
  
  numberGridCellsX = (int) ceil(range_width.x() / grid_cell_width + 0.5);
  numberGridCellsY = (int) ceil(range_width.y() / grid_cell_width + 0.5);
  numberGridCellsZ = (int) ceil(range_width.z() / grid_cell_width + 0.5);
  numberGridCells  = numberGridCellsX * numberGridCellsY * numberGridCellsZ;
  
  gridCellWidth.x(range_width.x() / ((double) numberGridCellsX));
  gridCellWidth.y(range_width.y() / ((double) numberGridCellsY));
  gridCellWidth.z(range_width.z() / ((double) numberGridCellsZ));
  
  gridCellWidthInverse.x(1.0 / gridCellWidth.x());
  gridCellWidthInverse.y(1.0 / gridCellWidth.y());
  gridCellWidthInverse.z(1.0 / gridCellWidth.z());
  
    // allocate neighborhood list array
  neighborhoodList = new DLIList<CubitPoint*>* [numberGridCells];
  assert(neighborhoodList != NULL);
  for ( i = 0; i < numberGridCells; i++)
  {
    neighborhoodList[i] = NULL;
  }
}

//=========================================================================== 
//Function Name: evaluate_position 
// 
//Member Type:  PUBLIC 
//Description:  evaluate the facet edge at a position 
//              eval_tangent is NULL if tangent not needed 
//=========================================================================== 
CubitStatus CubitFacetEdge::evaluate_position( const CubitVector &start_position, 
                                               CubitVector *eval_point, 
                                               CubitVector *eval_tangent)  
{ 
  CubitStatus stat = CUBIT_SUCCESS; 
 
  // find the adjacent facet 
 
  CubitFacet *facet_ptr = this->adj_facet( 0 ); 
 
  // If there is none or this is a linear representation -  
  // then project to the linear edge  
 
  if (!facet_ptr || facet_ptr->eval_order() == 0 || facet_ptr->is_flat()) 
  { 
    if (eval_point) 
    { 
      closest_point(start_position, *eval_point); 
    } 
    if (eval_tangent) 
    { 
      *eval_tangent = point(1)->coordinates() - point(0)->coordinates(); 
      (*eval_tangent).normalize(); 
    } 
  } 
  else 
  { 
    int vert0 = facet_ptr->point_index( point(0) ); 
    int vert1 = facet_ptr->point_index( point(1) ); 
    CubitVector pt_on_plane, close_point; 
    CubitVector start = start_position; 
    double dist_to_plane; 
    CubitBoolean outside_facet; 
    FacetEvalTool::project_to_facet_plane( facet_ptr, start,  
                                           pt_on_plane, dist_to_plane ); 
    stat = FacetEvalTool::project_to_facetedge( facet_ptr,  
                                                vert0, vert1, 
                                                start, 
                                                pt_on_plane,  
                                                close_point, 
                                                outside_facet ); 
    if (eval_point) 
    { 
      *eval_point = close_point; 
    } 
    if (eval_tangent) 
    { 
      CubitVector edvec = point(1)->coordinates() - point(0)->coordinates(); 
      edvec.normalize(); 
      CubitVector areacoord; 
      FacetEvalTool::facet_area_coordinate( facet_ptr, close_point, areacoord );  
      FacetEvalTool::eval_facet_normal(facet_ptr, areacoord, *eval_tangent); 
      CubitVector cross = edvec * *eval_tangent; 
      *eval_tangent = *eval_tangent * cross; 
      (*eval_tangent).normalize(); 
    } 
  } 
  return CUBIT_SUCCESS; 
} 

double PST_Edge::closest_on_line( const CubitVector& P )
{
  CubitVector B = start_coord();
  CubitVector M   = end_coord() - B;
  if( M.length_squared() < RESABS_SQR )
    return 0.0;
  return ( M % ( P - B ) ) / ( M % M );
}

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 );
} 

//-------------------------------------------------------------------------
// 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;
}

void PointGridSearch::get_neighborhood_points_sorted(
  DLIList<CubitPoint*> &point_list,
  const CubitVector& center, double cut_off)
{
  point_list.clean_out();
  DLIList<CubitPoint*> temp_point_list;
  int i;

  for (int k = boundingCellMinimumZ; k <= boundingCellMaximumZ; k++)
  {
    int kn = numberGridCellsY * k;
    for (int j = boundingCellMinimumY; j <= boundingCellMaximumY; j++)
    {
      int jn = numberGridCellsX * (kn + j);
      for ( i = boundingCellMinimumX; i <= boundingCellMaximumX; i++)
      {
	int in = jn + i;
        if (neighborhoodList[in])
        {
          temp_point_list += *(neighborhoodList[in]);
        }
      }
    }
  }

  // evaluate point distance to center ... remove those larger than cut_off
  SDLByDouble sorted_index_list;
  IndexedDouble *ID;  
  CubitVector vec;
  cut_off *= cut_off;
  temp_point_list.reset();
  for ( i = 0; i < temp_point_list.size(); i++)
  {
    vec = center - temp_point_list.get_and_step()->coordinates();
    double distance = vec.length_squared();
    if (distance < cut_off)
    {
      ID = new IndexedDouble( i, distance );
      sorted_index_list.append( ID );
    }
  }

  sorted_index_list.sort();
  temp_point_list.reset();
  for ( i = 0; i < sorted_index_list.size(); i++ )
  {
    ID = sorted_index_list.get_and_step();
    point_list.append( temp_point_list.next( ID->index() ) );
    delete ID;
  }
}

//===================================================================================
// 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       : 
//
// Special Notes : 
//
// Creator       : Jason Kraftcheck
//
// Creation Date : 05/10/04
//-------------------------------------------------------------------------
CubitStatus PartitionBody::mass_properties( CubitVector& result, double& volume )
{
  DLIList<Lump*> lump_list;
  lumps( lump_list );
  
  DLIList<PartitionLump*> part_list;
  CAST_LIST( lump_list, part_list, PartitionLump );
  if (part_list.size() < lump_list.size())
    return real_body()->mass_properties( result, volume );
  
  CubitVector centroid(0.0, 0.0, 0.0), tmp_centroid;
  volume = 0.0;
  double tmp_volume;
  for (int i = part_list.size(); i--; )
  {
    if (CUBIT_FAILURE == 
        part_list.get_and_step()->mass_properties( tmp_centroid, tmp_volume ))
      return CUBIT_FAILURE;
    
    centroid += tmp_volume * tmp_centroid;
    volume += tmp_volume;
  }
  
  if (volume > CUBIT_RESABS)
  {
    result = centroid / volume;
  }
  else
  {
    result.set( 0.0, 0.0, 0.0 );
    volume = 0.0;
  }
  return CUBIT_SUCCESS;
}

//===================================================================================
// Function: transform_to_uv (Public)
// Description: same as title, the local sizing will be returned in the z coord 
// Author: chynes
// Date: 7/10/02
//===================================================================================
CubitStatus PeriodicParamTool::transform_to_uv(const CubitVector &xyz_location, CubitVector &uv_location) 
{
	 
	double u,v;

	CubitStatus rv = refSurf->u_v_from_position(xyz_location, u, v);

	// offset values to avoid parameter discontinuity

	if (uPeriod && u < uOffset)
	{
		u += uPeriod;
	}
  
  // mirror surface if required to get correct loop orientation
  if (mirrorSurface)
  {
    u = -u;
  }

	if (vPeriod && v < vOffset)
	{
		v += vPeriod;
	}

	uv_location.set(u,v,1.0);


	return rv;
}