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

CubitStatus RefEdge::get_point_direction( CubitVector& origin, CubitVector& direction )
{
   Curve* curve_ptr = get_curve_ptr();

   if( curve_ptr != NULL )
   {

      if( curve_ptr->geometry_type() != STRAIGHT_CURVE_TYPE ) 
         return CUBIT_FAILURE;

      if( curve_ptr->get_point_direction( origin, direction ) == CUBIT_FAILURE )
      {
         origin = start_coordinates();
         direction = end_coordinates() - origin;
         direction.normalize();
      }
   }
   else 
   {
      PRINT_WARNING("In RefEdge::get_point_direction\n"
                    "         %s (curve %d) is not associated with a valid\n"
                    "         underlying geometric Curve\n",
                    entity_name().c_str(), id()) ;
      return CUBIT_FAILURE ;
   }
   
   if (curve_ptr->bridge_sense() == CUBIT_REVERSED)
     direction = -direction;
   
   return CUBIT_SUCCESS;

}

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

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

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

//====================================================================== 
// Function: edge_tangent 
// Description: return tangent at a point on the edge 
// Author: sjowen 
// Date: 2/01 
//====================================================================== 
CubitStatus  CubitFacetEdge::edge_tangent(  
  const CubitVector &/*point_on_edge*/,  
  CubitVector &tangent ) 
{ 
  CubitStatus stat = CUBIT_SUCCESS; 
  tangent = point(1)->coordinates() - 
            point(0)->coordinates(); 
  tangent.normalize(); 
  return stat; 
}

double CubitFacetEdge::angle_between_facets()
{
  CubitFacet *facet0 = adj_facet(0);
  CubitFacet *facet1 = adj_facet(1);

  // assumes facets are always oriented with outwards pointing normal

  CubitVector n0 = facet0->normal();
  CubitVector n1 = facet1->normal();

  // get orientation of edge with respect to facet0

  CubitPoint *p0 = point(0);
  CubitPoint *p1 = point(1);
  CubitPoint *pnext = facet0->next_node(p0);
  if (pnext != p1)
  {
    pnext = p0;
    p0 = p1;
    p1 = pnext;
  }
  CubitVector evec = p1->coordinates() - p0->coordinates();
  evec.normalize();
  CubitVector cross = n0 * n1;  cross.normalize();
  double angle;
  double edot = evec % cross;
  double ndot = n0 % n1;
  if (ndot >= 1.0)
    angle = 0.0;
  else if (ndot <= -1.0)
    angle = CUBIT_PI;
  else 
    angle = acos(ndot);
  if (edot <= 0.0)
  {
    angle = 2.0 * CUBIT_PI - angle;
  }
  return angle;
}

//=========================================================================== 
//Function Name: evaluate 
// 
//Member Type:  PUBLIC 
//Description:  evaluate the facet at area coordinates 
//              eval_normal is NULL if normal not needed 
//Note:         t is a value from -1 to 1.  t=0 is the midpoint 
//=========================================================================== 
CubitStatus CubitFacetEdge::evaluate( double &t,  
                                      CubitVector *eval_point, 
                                      CubitVector *eval_tangent )      
{ 
  CubitStatus stat = CUBIT_SUCCESS; 
 
  // project the position to the linear edge 
   
  double tt = (t + 1) * 0.5; 
  if (tt <= 0.0) tt = 0.0; 
  if (tt >= 1.0) tt = 1.0; 
  *eval_point = point(0)->coordinates() +  
    tt * (point(1)->coordinates() - point(0)->coordinates()); 
 
  // evaluate the point on the facet (if the order is higher than 0) 
 
  CubitFacet *facet_ptr = this->adj_facet( 0 ); 
  if (!facet_ptr || facet_ptr->is_flat()) 
  { 
    if (eval_tangent) 
    { 
      *eval_tangent = point(1)->coordinates() - point(0)->coordinates(); 
      (*eval_tangent).normalize(); 
    } 
  } 
  else 
  { 
    CubitVector areacoord; 
    FacetEvalTool::facet_area_coordinate( facet_ptr, *eval_point, areacoord ); 
    stat = facet_ptr->evaluate( areacoord, eval_point, eval_tangent ); 
    if (stat != CUBIT_SUCCESS) 
      return stat; 
    if (eval_tangent) 
    { 
      CubitVector edvec = point(1)->coordinates() - point(0)->coordinates(); 
      edvec.normalize(); 
      CubitVector cross = edvec * *eval_tangent; 
      *eval_tangent = *eval_tangent * cross; 
      (*eval_tangent).normalize(); 
    } 
  } 
  return stat; 
} 

//=============================================================================
//Function:  feature_angle (PRIVATE)
//Description: compute angles at nodes on the curve to see if we need to split
//             the curve.  Mark the node tooldata hitflag if the node will
//             break the curve (this is refernced in next_edge)
//Author: sjowen
//Date: 12/4/00
//=============================================================================
CubitStatus ChollaCurve::feature_angle(
  double min_dot )
{
  // first compute all of the edge vector and store with the edge tooldata

  int ii, jj;
  FacetEntity *facet_ptr;
  CubitFacetEdge *edge_ptr;
  CubitPoint *start_node;
  CubitPoint *end_node;
  CubitVector tangent;
  TDGeomFacet *td_gm;
  for (ii=0; ii<curveEdgeList.size(); ii++)
  {

    // compute the tangent vector of the edge and store it with its tooldata

    facet_ptr = curveEdgeList.get_and_step();
    edge_ptr = CAST_TO( facet_ptr, CubitFacetEdge );
    start_node = edge_ptr->point(0);
    end_node = edge_ptr->point(1);
    tangent = end_node->coordinates() - start_node->coordinates();
    tangent.normalize();
    td_gm = TDGeomFacet::get_geom_facet( edge_ptr );
    td_gm->set_normal( tangent );

    // initialize the nodes tooldata hit flags - set them all to -1

    td_gm = TDGeomFacet::get_geom_facet(start_node);
    td_gm->set_hit_flag(-1);
    td_gm = TDGeomFacet::get_geom_facet(end_node);
    td_gm->set_hit_flag(-1);
  }

  // now go through them again and compute the dot product between edges

  CubitVector tang0;
  CubitVector tang1;
  double dot;
  CubitPoint *node_ptr;
  CubitFacetEdge *next_edge_ptr;
  TDGeomFacet *td_gm_node;

  for (ii=0; ii<curveEdgeList.size(); ii++)
  {
    facet_ptr = curveEdgeList.get_and_step();
    edge_ptr = CAST_TO( facet_ptr, CubitFacetEdge );
    start_node = edge_ptr->point(0);
    end_node = edge_ptr->point(1);
    for (jj=0; jj<2; jj++)
    {
      node_ptr = (jj==0) ? start_node : end_node;
      td_gm_node = TDGeomFacet::get_geom_facet( node_ptr );
      if (td_gm_node->get_hit_flag() == -1)
      {
        next_edge_ptr = next_edge( node_ptr, edge_ptr );
        if (next_edge_ptr == NULL)
        {
          td_gm_node->set_hit_flag( 1 );
          node_ptr->set_as_feature();
        }
        else
        {
          td_gm = TDGeomFacet::get_geom_facet( edge_ptr );
          tang0 = td_gm->get_normal();
          td_gm = TDGeomFacet::get_geom_facet( next_edge_ptr );
          tang1 = td_gm->get_normal();

          // change the sign of the tangent vectors if the
          // sense of the edges are not the same

          if (node_ptr == start_node)
          {
            if (node_ptr != next_edge_ptr->point(1))
              tang0 = -tang0;
          }
          else
          {
            if (node_ptr != next_edge_ptr->point(0))
              tang0 = -tang0;
          }

          // compute the dot product between tangemt vectors

          dot = tang0 % tang1;

          // set the hit flag if there needs to be a feature break here

          if (dot <= min_dot)
          {
            td_gm_node->set_hit_flag( 1 );
            node_ptr->set_as_feature();
//.........这里部分代码省略.........

CubitStatus
SimplifyTool::weighted_average_normal(RefFace* ref_face,
                                      CubitVector &normal, 
                                      double &weight )
{
    GMem g_mem;
    unsigned short norm_tol = 30;
    double dist_tol = -1.0;

    ref_face->get_geometry_query_engine()->
        get_graphics(ref_face->get_surface_ptr(), &g_mem, norm_tol, dist_tol );

    if(g_mem.fListCount < 1)
    {
        // Decrease tolerance and try again (we can get this for small features)
        norm_tol /= 2;
        ref_face->get_geometry_query_engine()->
            get_graphics(ref_face->get_surface_ptr(), &g_mem, norm_tol, dist_tol );
    }

    if(g_mem.fListCount < 1)
    {
        // Lets give up 
        PRINT_ERROR( "Unable to find average normal of a surface\n" );
        return CUBIT_FAILURE;
    }

    // Initialize
    weight = 0.0;
    normal.set( 0.0, 0.0, 0.0 );

    // Loop through the triangles
    double tri_weight, A, B, C;
    GPoint p[3];
    GPoint* plist = g_mem.point_list();
    int* facet_list = g_mem.facet_list();
    int c = 0;
    for( ;c<g_mem.fListCount; )
    {
        p[0] = plist[facet_list[++c]];
        p[2] = plist[facet_list[++c]];
        p[1] = plist[facet_list[++c]]; 
        c++;

        // Get centroid
        CubitVector p1( p[0].x, p[0].y, p[0].z );
        CubitVector p2( p[2].x, p[2].y, p[2].z );
        CubitVector p3( p[1].x, p[1].y, p[1].z );

        CubitVector center = (p1 + p2 + p3)/3.0;

        CubitVector norm(ref_face->normal_at(center));

        // Get triangle area
        A = p1.y() * p2.z() + p1.z() * p3.y() + p2.y() * p3.z() -
            p2.z() * p3.y() - p1.y() * p3.z() - p1.z() * p2.y();

        B = p1.z() * p2.x() + p1.x() * p3.z() + p2.z() * p3.x() -
            p2.x() * p3.z() - p1.z() * p3.x() - p1.x() * p2.z();

        C = p1.x() * p2.y() + p1.y() * p3.x() + p2.x() * p3.y() -
            p2.y() * p3.x() - p1.x() * p3.y() - p1.y() * p2.x();

        //Note: triangle area = 0.5*(sqrt(A*A+B*B+C*C));

        tri_weight = 0.5*(A*A+B*B+C*C);

        normal += tri_weight * norm;

        weight += tri_weight;

    }

    normal.normalize();

    return CUBIT_SUCCESS;
}