C++ VerdictVector::length_squared方法代码示例

/*!
   the radius ratio of a triangle

   NB (P. Pebay 01/13/07): 
     CR / (2.0*IR) where CR is the circumradius and IR is the inradius

     The radius ratio is also known to VERDICT, for tetrahedral elements only,
     as the "aspect beta".
   
*/
C_FUNC_DEF double v_tri_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
{

  // three vectors for each side 
  VerdictVector a( coordinates[1][0] - coordinates[0][0],
                   coordinates[1][1] - coordinates[0][1],
                   coordinates[1][2] - coordinates[0][2] );
  
  VerdictVector b( coordinates[2][0] - coordinates[1][0],
                   coordinates[2][1] - coordinates[1][1],
                   coordinates[2][2] - coordinates[1][2] );
  
  VerdictVector c( coordinates[0][0] - coordinates[2][0],
                   coordinates[0][1] - coordinates[2][1],
                   coordinates[0][2] - coordinates[2][2] );

  double a1 = a.length();
  double b1 = b.length();
  double c1 = c.length();
 
  VerdictVector ab = a * b;
  double denominator = ab.length_squared();

  if( denominator < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;

  double radius_ratio;
  radius_ratio = .25 * a1 * b1 * c1 * ( a1 + b1 + c1 ) / denominator;
  
  if( radius_ratio > 0 )
    return (double) VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
}

inline double normalize_jacobian( double jacobi,
                                  VerdictVector& v1,
                                  VerdictVector& v2,
                                  VerdictVector& v3,
                                  int tet_flag = 0 )
{
    double return_value = 0.0;

    if ( jacobi != 0.0 )
    {

        double l1, l2, l3, length_product;
        // Note: there may be numerical problems if one is a lot shorter
        // than the others this way. But scaling each vector before the
        // triple product would involve 3 square roots instead of just
        // one.
        l1 = v1.length_squared();
        l2 = v2.length_squared();
        l3 = v3.length_squared();
        length_product = sqrt( l1 * l2 * l3 );

        // if some numerical scaling problem, or just plain roundoff,
        // then push back into range [-1,1].
        if ( length_product < fabs(jacobi) )
        {
            length_product = fabs(jacobi);
        }

        if( tet_flag == 1)
            return_value = v_sqrt_2 * jacobi / length_product;
        else
            return_value = jacobi / length_product;

    }
    return return_value;

}

/*! 
  tri_quality for calculating multiple tri functions at once

  using this method is generally faster than using the individual 
  method multiple times.

*/
C_FUNC_DEF void v_tri_quality( int num_nodes, double coordinates[][3], 
    unsigned int metrics_request_flag, TriMetricVals *metric_vals ) 
{

  memset( metric_vals, 0, sizeof(TriMetricVals) );

  // for starts, lets set up some basic and common information

  /*  node numbers and side numbers used below

             2
             ++
            /  \ 
         2 /    \ 1
          /      \
         /        \
       0 ---------+ 1
             0
  */
  
  // vectors for each side
  VerdictVector sides[3];
  sides[0].set(
      coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2]
      );
  sides[1].set(
      coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2]
      );
  sides[2].set(
      coordinates[2][0] - coordinates[0][0],
      coordinates[2][1] - coordinates[0][1],
      coordinates[2][2] - coordinates[0][2]
      );
  VerdictVector tri_normal = sides[0] * sides[2];
    //if we have access to normal information, check to see if the
    //element is inverted.  If we don't have the normal information
    //that we need for this, assume the element is not inverted.
    //This flag will be used for condition number, jacobian, shape,
    //and size and shape.
  bool is_inverted = false;
  if( compute_normal )
  {
      //center of tri
    double point[3], surf_normal[3];
    point[0] =  (coordinates[0][0] + coordinates[1][0] + coordinates[2][0]) / 3;
    point[1] =  (coordinates[0][1] + coordinates[1][1] + coordinates[2][1]) / 3;
    point[2] =  (coordinates[0][2] + coordinates[1][2] + coordinates[2][2]) / 3;
      //dot product
    compute_normal( point, surf_normal ); 
    if( (tri_normal.x()*surf_normal[0] + 
         tri_normal.y()*surf_normal[1] +
         tri_normal.z()*surf_normal[2] ) < 0 )
      is_inverted=true; 
  }
  
  // lengths squared of each side
  double sides_lengths_squared[3];
  sides_lengths_squared[0] = sides[0].length_squared();
  sides_lengths_squared[1] = sides[1].length_squared();
  sides_lengths_squared[2] = sides[2].length_squared();
 

  // if we are doing angle calcuations
  if( metrics_request_flag & (V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE) )
  {
    // which is short and long side
    int short_side=0, long_side=0;

    if(sides_lengths_squared[1] < sides_lengths_squared[0])
      short_side = 1;
    if(sides_lengths_squared[2] < sides_lengths_squared[short_side])
      short_side = 2;
  
    if(sides_lengths_squared[1] > sides_lengths_squared[0])
      long_side = 1;
    if(sides_lengths_squared[2] > sides_lengths_squared[long_side])
      long_side = 2;


    // calculate the minimum angle of the tri
    if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
    {
      if(sides_lengths_squared[0] == 0.0 || 
        sides_lengths_squared[1] == 0.0 ||
        sides_lengths_squared[2] == 0.0)
      {
        metric_vals->minimum_angle = 0.0;
      }        
      else if(short_side == 0)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      scaled_jac = areas[1] / (lengths[1] * lengths[0]);
      min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

      scaled_jac = areas[2] / (lengths[2] * lengths[1]);
      min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

      scaled_jac = areas[3] / (lengths[3] * lengths[2]);
      min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

      metric_vals->scaled_jacobian = min_scaled_jac;
     
      //what the heck...set shear as well 
      if( min_scaled_jac <= VERDICT_DBL_MIN )
        metric_vals->shear = 0.0;
      else
        metric_vals->shear = min_scaled_jac;
    }
  }

  if( metrics_request_flag & (V_QUAD_WARPAGE | V_QUAD_ODDY) )
  {
    VerdictVector corner_normals[4];
    corner_normals[0] = edges[3] * edges[0];
    corner_normals[1] = edges[0] * edges[1];
    corner_normals[2] = edges[1] * edges[2];
    corner_normals[3] = edges[2] * edges[3];

    if( metrics_request_flag & V_QUAD_ODDY )
    {
      double oddy, max_oddy = 0.0;

      double diff, dot_prod;

      double length_squared[4];
      length_squared[0] = corner_normals[0].length_squared();
      length_squared[1] = corner_normals[1].length_squared();
      length_squared[2] = corner_normals[2].length_squared();
      length_squared[3] = corner_normals[3].length_squared();

      if( length_squared[0] < VERDICT_DBL_MIN ||
          length_squared[1] < VERDICT_DBL_MIN ||
          length_squared[2] < VERDICT_DBL_MIN ||
          length_squared[3] < VERDICT_DBL_MIN )
       metric_vals->oddy = VERDICT_DBL_MAX; 
      else
      {
        diff = (lengths[0]*lengths[0]) - (lengths[1]*lengths[1]);
        dot_prod = edges[0]%edges[1];
        oddy = ((diff*diff) + 4*dot_prod*dot_prod ) / (2*length_squared[1]);
        max_oddy = VERDICT_MAX( oddy, max_oddy );

        diff = (lengths[1]*lengths[1]) - (lengths[2]*lengths[2]);
        dot_prod = edges[1]%edges[2];
        oddy = ((diff*diff) + 4*dot_prod*dot_prod ) / (2*length_squared[2]);
        max_oddy = VERDICT_MAX( oddy, max_oddy );

        diff = (lengths[2]*lengths[2]) - (lengths[3]*lengths[3]);
        dot_prod = edges[2]%edges[3];
        oddy = ((diff*diff) + 4*dot_prod*dot_prod ) / (2*length_squared[3]);
        max_oddy = VERDICT_MAX( oddy, max_oddy );

        diff = (lengths[3]*lengths[3]) - (lengths[0]*lengths[0]);
        dot_prod = edges[3]%edges[0];
        oddy = ((diff*diff) + 4*dot_prod*dot_prod ) / (2*length_squared[0]);
        max_oddy = VERDICT_MAX( oddy, max_oddy );