cgma
LoopParamTool.cpp
Go to the documentation of this file.
00001 //- Class: LoopParamTool
00002 //-------------------------------------------------------------------------
00003 // Filename      : LoopParamTool.cpp
00004 //
00005 // Purpose       : Surface Parameterization for mesh by triangulation flattening
00006 //
00007 // Creator       : Shiraj Khan
00008 //
00009 // Creation Date : 1/21/2003
00010 //
00011 // Owner         : Shiraj Khan
00012 //-------------------------------------------------------------------------
00013 
00014 #include "LoopParamTool.hpp"
00015 #include "CastTo.hpp"
00016 #include "CubitPointData.hpp"
00017 #include "GeometryDefines.h"
00018 #include "GeometryDefines.h"
00019 #include "CubitMessage.hpp"
00020 #include "math.h"
00021 #define EPSILON_LOWER -0.00000001
00022 #define EPSILON_UPPER  0.00000001
00023 #define SWAP(a,b) { temp = a; a = b; b = temp; }
00024 
00025 
00026 
00027 
00028 //-------------------------------------------------------------------------
00029 // Function:    LoopParamTool
00030 // Description: constructor
00031 // Author:      Shiraj Khan
00032 // Date:        1/21/2003
00033 //-------------------------------------------------------------------------
00034 LoopParamTool::LoopParamTool()
00035 {
00036     
00037 }
00038 
00039 //-------------------------------------------------------------------------
00040 // Function:    LoopParamTool
00041 // Description: deconstructor
00042 // Author:      Shiraj Khan
00043 // Date:        1/21/2003
00044 //-------------------------------------------------------------------------
00045 LoopParamTool::~LoopParamTool() {}
00046 
00047 //===================================================================================
00048 // Function: set_up_space (Public)
00049 // Description: sets up space of flattening
00050 // Author: Shiraj Khan
00051 // Date: 1/21/2003
00052 //===================================================================================
00053 CubitStatus LoopParamTool::set_up_space() 
00054 {
00055   return CUBIT_SUCCESS;
00056 }
00057 
00058 //===================================================================================
00059 // Function: new_space_LoopParam (Public)
00060 // Description: sets up space of flattening
00061 // Author: Shiraj Khan
00062 // Date: 1/21/2003
00063 //===================================================================================
00064 CubitStatus LoopParamTool::new_space_LoopParam( DLIList<DLIList<CubitPoint *>*> &loops_cubit_points,
00065                                                 CubitVector* normal)
00066 {
00067     
00068   int ii;
00069   int jj;
00070   double sumx = 0.0, sumy = 0.0, sumz = 0.0;
00071   double sumxx = 0.0, sumxy = 0.0, sumxz = 0.0;
00072   double sumyy = 0.0, sumyz = 0, sumzz = 0.0;
00073   double aa[4][4];
00074   double Xc = 0.0, Yc = 0.0, Zc = 0.0;
00075   int n = 0;
00076   
00077   CubitPoint *point, *point1, *point2;
00078   CubitVector point_coordinates, point_coordinates1, point_coordinates2;
00079   CubitVector cm_plane; // center of mass of a plane
00080   CubitVector vec1, vec2, plane_normal;
00081   DLIList<CubitPoint *> *sub_loops_cubit_points;
00082 
00083   
00084   for ( ii = 0; ii < loops_cubit_points.size(); ii++ )
00085   {
00086     sub_loops_cubit_points = loops_cubit_points.get_and_step();
00087     for ( jj = 0; jj < sub_loops_cubit_points->size(); jj++ )
00088     {
00089       point = sub_loops_cubit_points->get_and_step();
00090       point_coordinates = point->coordinates();
00091       sumx = sumx + point_coordinates.x();
00092       sumy = sumy + point_coordinates.y();
00093       sumz = sumz + point_coordinates.z();
00094       sumxx = sumxx + ( point_coordinates.x() * point_coordinates.x() );
00095       sumxy = sumxy + ( point_coordinates.x() * point_coordinates.y() );
00096       sumxz = sumxz + ( point_coordinates.x() * point_coordinates.z() );
00097       sumyy = sumyy + ( point_coordinates.y() * point_coordinates.y() );
00098       sumyz = sumyz + ( point_coordinates.y() * point_coordinates.z() );
00099       sumzz = sumzz + ( point_coordinates.z() * point_coordinates.z() );
00100       n++;
00101     }
00102   }
00103   Xc = sumx / n;
00104   Yc = sumy / n;
00105   Zc = sumz / n;
00106   cm_plane.set(Xc,Yc,Zc);
00107   if ( Xc < EPSILON_UPPER && Xc > EPSILON_LOWER ) Xc = 0.0;
00108   if ( Yc < EPSILON_UPPER && Yc > EPSILON_LOWER ) Yc = 0.0;
00109   if ( Zc < EPSILON_UPPER && Zc > EPSILON_LOWER ) Zc = 0.0;
00110   aa[1][1] = sumxx - Xc * sumx;
00111   aa[1][2] = sumxy - Yc * sumx;
00112   aa[1][3] = sumxz - Zc * sumx;
00113   aa[2][1] = sumxy - Xc * sumy;
00114   aa[2][2] = sumyy - Yc * sumy;
00115   aa[2][3] = sumyz - Zc * sumy;
00116   aa[3][1] = sumxz - Xc * sumz;
00117   aa[3][2] = sumyz - Yc * sumz;
00118   aa[3][3] = sumzz - Zc * sumz;
00119     
00120   for ( ii = 1; ii <=3; ii++ )
00121   {
00122     for ( jj = 1; jj <=3; jj++ )
00123       if ( aa[ii][jj] < EPSILON_UPPER && aa[ii][jj] > EPSILON_LOWER )
00124         aa[ii][jj] = 0.0;
00125   }
00126   double determinant_aa = aa[1][1] * ( aa[2][2]*aa[3][3] - aa[3][2]*aa[2][3] ) - aa[1][2] * ( aa[2][1]*aa[3][3] - aa[3][1]*aa[2][3] )
00127     + aa[1][3] * ( aa[2][1]*aa[3][2] - aa[3][1]*aa[2][2] );
00128 
00129   if ( determinant_aa < EPSILON_UPPER && determinant_aa > EPSILON_LOWER ) 
00130     determinant_aa = 0.0;
00131   loops_cubit_points.reset();
00132   
00133     // if determinant is 0.0 ( all the points are lying on the plane), the equation of a plane: 
00134     //(vec1) crossproduct (vec2)
00135     // where vec1 and vec2 are the vectors originating from a common point( center of mass of a plane) and  
00136     // lying on the plane.
00137   if(normal){
00138     a = normal->x();
00139     b = normal->y();
00140     c = normal->z();
00141     d = -(a*Xc + b*Yc + c*Zc );
00142   }
00143   else if ( determinant_aa == 0.0 )
00144   {
00145     sub_loops_cubit_points = loops_cubit_points.get_and_step();
00146     point1 = sub_loops_cubit_points->get_and_step();
00147     point2 = sub_loops_cubit_points->get_and_step();
00148     point_coordinates1 = point1->coordinates(); 
00149     point_coordinates2 = point2->coordinates();
00150     vec1 = cm_plane - point_coordinates1;
00151     vec2 = cm_plane - point_coordinates2;
00152     plane_normal = vec1 * vec2;
00153     a = plane_normal.x();
00154     b = plane_normal.y();
00155     c = plane_normal.z();
00156     d = -( a*Xc + b*Yc + c*Zc );
00157   }
00158   else // to calculate eigen vector corresponding to the smallest eigen value
00159   {
00160       // to calculate the inverse of a matrix
00161     int i, j, k;
00162     int icol = -1, irow = -1, l, ll;
00163     double big, z, pivinv, temp;
00164     int indxc[4];
00165     int indxr[4];
00166     int ipiv[4];
00167     for ( j = 1; j <= 3; j++)
00168     {
00169       indxc[j] = 0;
00170       indxr[j] = 0;
00171       ipiv[j] = 0;
00172     }
00173     for ( i = 1; i <= 3;i++)
00174     {
00175       big = 0.0;
00176       for ( j = 1; j <= 3; j++)
00177       {
00178         if ( ipiv[j] != 1 )
00179         {
00180           for ( k = 1; k <= 3; k++)
00181           {
00182             if (ipiv[k] == 0 )
00183             {
00184               if (fabs(aa[j][k]) >= big )
00185               {
00186                 big = fabs(aa[j][k]);
00187                 irow = j;
00188                 icol = k;
00189               }
00190             }
00191             else if(ipiv[k]>1) PRINT_ERROR("Matrix is singular1\n");
00192           }
00193         }
00194       }
00195       ipiv[icol]+=1;
00196       if ( irow != icol )
00197       {
00198         for (l=1; l <= 3;l++)
00199           SWAP( aa[irow][l],aa[icol][l] );
00200             
00201       }
00202       indxr[i] = irow;
00203       indxc[i] = icol;
00204       if(aa[icol][icol] == 0.0 ) PRINT_ERROR("Matrix is singular2\n");
00205       pivinv = 1.0/aa[icol][icol];
00206       aa[icol][icol] = 1.0;
00207       for( l = 1; l <= 3; l++)
00208         aa[icol][l] *= pivinv;
00209       for ( ll = 1; ll <= 3; ll++)
00210       {
00211         if ( ll != icol )
00212         {
00213           z = aa[ll][icol];
00214           aa[ll][icol] = 0.0;
00215           for ( l = 1;l <= 3; l++)
00216             aa[ll][l] -= aa[icol][l] * z;
00217                 
00218         }
00219       }
00220             
00221       for ( l = 3; l >= 1; l-- )
00222       {
00223         if ( indxr[l] != indxc[l] )
00224           for ( k = 1; k <= 3; k++ )
00225             SWAP( aa[k][indxr[l]], aa[k][indxc[l]] )
00226       }
00227     }
00228         
00229         
00230       // Power Method to get the Eigen Vector corresponding to the smallest Eigen value
00231     double x[4] = {0.0, 1.0, 1.1, 1.2};
00232     double dd[4];
00233     double zz[4];
00234     int t = 0;
00235     double diff1 = 0.0, diff2 = 0.0;
00236     int flag;
00237     double largest_dd;
00238     while ( t <= 100 )
00239     {
00240       flag = 1;
00241       for ( i = 1; i <= 3; i++ )
00242       {
00243         dd[i] = 0.0;
00244         for ( j = 1; j <= 3; j++ )
00245         {
00246           dd[i] = dd[i] + aa[i][j]  * x[j];
00247 
00248           if (dd[i] > EPSILON_LOWER && dd[i] < EPSILON_UPPER) dd[i] = 0.0;
00249         }
00250       }
00251       t = t+1;
00252       largest_dd = dd[1] > dd[2] ? ( dd[1] > dd[3] ? dd[1] : dd[3] ) : ( dd[2] > dd[3] ? dd[2] : dd[3] );
00253       for ( i = 1; i <= 3; i++ )
00254       {
00255         zz[i] = dd[i]/largest_dd;
00256         if (zz[i] > EPSILON_LOWER && zz[i] < EPSILON_UPPER) zz[i] = 0.0;
00257       }
00258       for ( i = 1; i <= 3; i++ )
00259       {
00260         diff1 =  fabs(x[i] - zz[i]);
00261         diff2 = fabs(diff1) - EPSILON_UPPER*fabs(zz[i]);
00262         if ( diff2 <= 0.0 )
00263           continue;
00264         else
00265         {
00266           flag = 0;
00267           break;
00268         }
00269       }
00270       if ( flag == 1 )
00271         break;
00272       if (flag == 0 )
00273       {
00274         for ( i = 1; i <= 3; i++ )
00275           x[i] = zz[i];
00276       }
00277         
00278 
00279     }
00280     for ( i = 1; i <= 3; i++ )
00281     {
00282       x[i] = zz[i];
00283       if ( (x[i] < EPSILON_UPPER) && (x[i] > EPSILON_LOWER) )
00284         x[i] = 0.0;
00285     }
00286     a = x[1];
00287     b = x[2];
00288     c = x[3];
00289     d = -(a*Xc+b*Yc+c*Zc);
00290         
00291   }
00292 
00293   //printf("\n\nEquation of Plane is %fx+%fy+%fz+%f = 0\n",a,b,c,d);
00294   loops_cubit_points.reset();
00295   sub_loops_cubit_points = loops_cubit_points.get_and_step();
00296   CubitPoint *p1 = sub_loops_cubit_points->get_and_step();
00297   CubitVector p1_coordinates = p1->coordinates();
00298   double p1x = p1_coordinates.x() - (p1_coordinates.x()*a + p1_coordinates.y()*b + p1_coordinates.z()*c + d)*
00299     (a)/(pow(a,2)+pow(b,2)+pow(c,2));
00300   double p1y = p1_coordinates.y() - (p1_coordinates.x()*a + p1_coordinates.y()*b + p1_coordinates.z()*c + d)*
00301     (b)/(pow(a,2)+pow(b,2)+pow(c,2));
00302   double p1z = p1_coordinates.z() - (p1_coordinates.x()*a + p1_coordinates.y()*b + p1_coordinates.z()*c + d)*
00303     (c)/(pow(a,2)+pow(b,2)+pow(c,2));
00304     
00305   p1_coordinates.set(p1x, p1y, p1z);
00306   CubitVector surf_normal;
00307   surf_normal.x(a); surf_normal.y(b); surf_normal.z(-1.0);
00308   surf_normal.normalize();
00309 
00310     // use the first two points on the boundary to define gradient vectors
00311     // and orient our uv space.
00312   CubitPoint *p2;
00313   CubitVector p2_coordinates;
00314   double p2x, p2y, p2z;
00315   double cos = 1.0;
00316     // 
00317   while (cos == 1.0)
00318   {
00319     p2 = sub_loops_cubit_points->get_and_step();
00320     p2_coordinates = p2->coordinates();
00321     p2x = p2_coordinates.x() - (p2_coordinates.x()*a + p2_coordinates.y()*b + p2_coordinates.z()*c + d)*
00322       (a)/(pow(a,2)+pow(b,2)+pow(c,2));
00323     p2y = p2_coordinates.y() - (p2_coordinates.x()*a + p2_coordinates.y()*b + p2_coordinates.z()*c + d)*
00324       (b)/(pow(a,2)+pow(b,2)+pow(c,2));
00325     p2z = p2_coordinates.z() - (p2_coordinates.x()*a + p2_coordinates.y()*b + p2_coordinates.z()*c + d)*
00326       (c)/(pow(a,2)+pow(b,2)+pow(c,2));
00327     p2_coordinates.set(p2x, p2y, p2z);
00328         
00329         
00330     if ( p1_coordinates == p2_coordinates )
00331       cos = 1.0;
00332     else cos = 0.0;
00333   }
00334   Du = p2_coordinates - p1_coordinates;
00335   Du.normalize();
00336   Dv = surf_normal * Du;
00337   Dv.normalize();
00338   uvCenter = p1_coordinates;
00339   sub_loops_cubit_points->reset();
00340 
00341   return CUBIT_SUCCESS; 
00342 }
00343 
00344 //===================================================================================
00345 // Function: transform_to_bestfit_plane (Public)
00346 // Description: same as title
00347 // Author: Shiraj Khan
00348 // Date: 1/21/2003
00349 //===================================================================================
00350 CubitStatus LoopParamTool::transform_to_bestfit_plane(DLIList<DLIList<CubitPoint *>*> &loop_cubit_points)
00351 {
00352   int ii, jj;
00353   DLIList<CubitPoint *> *sub_loop_cubit_points = NULL;
00354   CubitPoint *point_ptr = NULL;
00355     
00356   CubitVector point_coordinates;
00357   double x, y, z;
00358 
00359   for ( ii = 0; ii < loop_cubit_points.size(); ii++ )
00360   {
00361     sub_loop_cubit_points = loop_cubit_points.get_and_step();
00362         
00363     for ( jj = 0; jj < sub_loop_cubit_points->size(); jj++ )
00364     {
00365       point_ptr = sub_loop_cubit_points->get_and_step();
00366       ToolData* td = point_ptr->get_TD( &TDVector::is_td_vector );
00367       TDVector* td_pos = NULL;
00368       if( NULL == td )
00369           td_pos = new TDVector(point_ptr->coordinates());
00370       else
00371           td_pos = static_cast<TDVector*>( td );
00372       point_ptr->add_TD(td_pos);
00373       point_coordinates = point_ptr->coordinates();
00374             
00375       x = point_coordinates.x() - (point_coordinates.x()*a + point_coordinates.y()*b + point_coordinates.z()*c + d)*
00376         (a)/(pow(a,2)+pow(b,2)+pow(c,2));
00377       y = point_coordinates.y() - (point_coordinates.x()*a + point_coordinates.y()*b + point_coordinates.z()*c + d)*
00378         (b)/(pow(a,2)+pow(b,2)+pow(c,2));
00379       z = point_coordinates.z() - (point_coordinates.x()*a + point_coordinates.y()*b + point_coordinates.z()*c + d)*
00380         (c)/(pow(a,2)+pow(b,2)+pow(c,2));
00381       point_coordinates.set(x, y, z);
00382       point_ptr->set(point_coordinates);
00383             
00384             
00385     }
00386         
00387   }
00388     
00389   return CUBIT_SUCCESS;
00390 }
00391 //===================================================================================
00392 // Function: transform_to_uv (Public)
00393 // Description: Does nothing for this tool.
00394 // Author: Shiraj Khan
00395 // Date: 1/21/2003
00396 //===================================================================================
00397 CubitStatus LoopParamTool::transform_to_uv(const CubitVector &, CubitVector &) 
00398 {
00399   PRINT_ERROR("This function is not appropriate for the the LoopParamTool.\n");
00400   assert(0);
00401   return CUBIT_FAILURE;
00402 }
00403 
00404 
00405 //===================================================================================
00406 // Function: transform_to_uv_local (Private)
00407 // Description: It transforms points from the Best Fit plane plane to UV space
00408 // Author: Shiraj Khan
00409 // Date: 1/21/2003
00410 //===================================================================================
00411 CubitStatus LoopParamTool::transform_to_uv_local(CubitVector &xyz_location, CubitVector &uv_location) 
00412 {
00413     // Translate to local origin at center
00414 
00415   CubitVector vect = xyz_location - uvCenter;
00416 
00417     // Multiply by transpose (inverse) of transformation vector
00418 
00419   uv_location.x( vect % Du );
00420   uv_location.y( vect % Dv );
00421   uv_location.z( 1.0 );
00422 
00423   return CUBIT_SUCCESS;
00424 }
00425 //===================================================================================
00426 // Function: transform_to_xyz (Private)
00427 // Description: Does nothing, just fills in virtual function.
00428 // Author: Shiraj Khan
00429 // Date: 1/21/2003
00430 //===================================================================================
00431 CubitStatus LoopParamTool::transform_to_xyz(CubitVector &, const CubitVector &) 
00432 {
00433   PRINT_ERROR("This function is not appropriate for the the LoopParamTool.\n");
00434   assert(0);
00435   return CUBIT_FAILURE;
00436 }
00437 
00438 //===================================================================================
00439 // Function: transform_to_xyz_local (Priavate)
00440 // Description: same as title
00441 // Author: Shiraj Khan
00442 // Date: 1/21/2003
00443 //===================================================================================
00444 CubitStatus LoopParamTool::transform_to_xyz_local(CubitVector &xyz_location, CubitVector &uv_location) 
00445 {
00446 // Multiply by transformation matrix
00447 
00448   CubitVector vect;
00449   vect.x( uv_location.x() * Du.x() +
00450           uv_location.y() * Dv.x() );
00451   vect.y( uv_location.x() * Du.y() +
00452           uv_location.y() * Dv.y() );
00453   vect.z( uv_location.x() * Du.z() +
00454           uv_location.y() * Dv.z() );
00455 
00456     // Translate from origin
00457 
00458   xyz_location = vect + uvCenter;
00459 
00460   return CUBIT_SUCCESS;
00461 }
00462 
00463 
00464 //-------------------------------------------------------------------------
00465 // Function:    transform_loopspoints_to_uv
00466 // Description: transform all the boundary points to the UV space.  Use
00467 //              the x-y coordinates in the points and set the z coordinate
00468 //              to one
00469 // Author:      Shiraj Khan
00470 // Date:        1/22/2003
00471 //-------------------------------------------------------------------------
00472 
00473 CubitStatus LoopParamTool::transform_loopspoints_to_uv( DLIList<DLIList<CubitPoint *>*> &loops_point_list)
00474                                                          
00475 {
00476   int ii, jj;
00477   CubitStatus ret_value = CUBIT_SUCCESS;
00478   DLIList<CubitPoint *> *sub_loops_point_list;
00479   CubitPoint *point_ptr;
00480   CubitVector xyz_location, uv_location, point_coordinates;
00481   
00482     // transform the points on the best fit plane
00483   ret_value = transform_to_bestfit_plane( loops_point_list);
00484   
00485     /*printf("\n after transform to best fit plane\n");
00486       for( ii = 0; ii < loops_point_list.size(); ii++)
00487       {
00488       sub_loops_point_list = loops_point_list.get_and_step();
00489       for ( jj = 0; jj < sub_loops_point_list->size(); jj++)
00490       {
00491       point_ptr = sub_loops_point_list->get_and_step();
00492       point_coordinates = point_ptr->coordinates();
00493       printf(" %f %f %f \n",point_coordinates.x(),point_coordinates.y(),point_coordinates.z());
00494       }
00495       }*/
00496 
00497     //transform to u-v
00498   for ( ii = 0; ii < loops_point_list.size(); ii++)
00499   {
00500     sub_loops_point_list = loops_point_list.get_and_step();
00501       
00502     for ( jj = 0; jj < sub_loops_point_list->size(); jj++ )
00503     {
00504       point_ptr = sub_loops_point_list->get_and_step();
00505       xyz_location = point_ptr->coordinates();
00506       ret_value = transform_to_uv_local( xyz_location, uv_location );
00507       point_ptr->set(uv_location);
00508           
00509     } 
00510   }
00511 
00512   if ( check_selfintersecting_coincident_edges(loops_point_list) != CUBIT_SUCCESS )
00513   {
00514     //PRINT_ERROR("Self interesections in parameter space found.\n");
00515     ret_value = CUBIT_FAILURE;
00516   }
00517   else
00518     ret_value = CUBIT_SUCCESS;
00519   
00520  
00521     // if the self-intersecting and coincident edges are not there then return CUBIT_SUCCESS
00522     // else restore the points to their original positions and return CUBIT_FAILURE
00523   return ret_value;
00524     //don't do the rest for now...
00525   if ( ret_value == CUBIT_SUCCESS )
00526     return CUBIT_SUCCESS;
00527   else
00528   {
00529     ToolData *td;
00530     for ( ii = 0; ii < loops_point_list.size(); ii++ )
00531     {
00532       sub_loops_point_list = loops_point_list.next(ii);
00533       for ( jj = 0; jj < sub_loops_point_list->size(); jj++)
00534       {
00535         point_ptr = sub_loops_point_list->next(jj);
00536         td = point_ptr->get_TD(&TDVector::is_td_vector);
00537         TDVector *td_pos = CAST_TO(td, TDVector);
00538         CubitVector orig_xyz_location = td_pos->get_vector();
00539         point_ptr->set(orig_xyz_location);
00540         point_coordinates = point_ptr->coordinates();
00541       }
00542     }
00543     return CUBIT_FAILURE;
00544   }
00545 }
00546 //-------------------------------------------------------------------------
00547 // Function:    check for self-intersecting and coincident edges
00548 // Description: check for self-intersecting and coincident edges. If the condition is true,
00549 //               it generates an error " Can't mesh due to self intersecting/coincident edges
00550 // Author:      Shiraj Khan
00551 // Date:        02/22/2003
00552 //-------------------------------------------------------------------------
00553 
00554 CubitStatus LoopParamTool::check_selfintersecting_coincident_edges(DLIList<DLIList<CubitPoint *>*> 
00555                                                                    &loops_bounding_points)
00556 {
00557     //This function checks each line segment against every other line segment
00558     //to see if they intersect or overlap.  Each line will be described
00559     //by the following equations:
00560     //
00561     // Pa = P1 + ua( P2 - P1 )
00562     // Pb = P3 + ub( P4 - P3 )
00563     //
00564     //Setting the Pa equal to Pb and solving for ua and ub you will get two
00565     //equation.  Both will have the common denominator:
00566     //
00567     // denom = ( (uv4.y() - uv3.y()) * (uv2.x() - uv1.x() ) -
00568     //         (uv4.x() - uv3.x()) * (uv2.y() - uv1.y() ) )
00569     //
00570     //And the numerators will be:
00571     //
00572     // numer_a = ( (uv4.x() - uv3.x()) * (uv1.y() - uv3.y() ) -
00573     //          (uv4.y() - uv3.y()) * (uv1.x() - uv3.x() ) )
00574     //      
00575     // numer_b = ( (uv2.x() - uv1.x()) * (uv1.y() - uv3.y() ) -
00576     //           (uv2.y() - uv1.y()) * (uv1.x() - uv3.x() ) )
00577     //
00578     //ua and ub then become:
00579     //
00580     // ua = numer_a / denom
00581     // ub = numer_b / denom
00582     //
00583     //If the lines are parallel then denom will be zero.  If they are
00584     //also coincedent then the numerators will also be zero.  For the
00585     //segments to intersect ua and ub need to both be between 0
00586     //and 1.
00587   
00588   int ii, jj, kk, ll;
00589   
00590   for ( ii = 0; ii < loops_bounding_points.size(); ii++ )
00591   {
00592     DLIList<CubitPoint*>* sub_loops_bounding_points1 = loops_bounding_points.next(ii);
00593     for( jj = 0; jj < sub_loops_bounding_points1->size(); jj++ )
00594     {
00595       CubitPoint* start_point1 = sub_loops_bounding_points1->next(jj);
00596       CubitPoint* end_point1 = NULL;
00597       if ( jj == sub_loops_bounding_points1->size()-1 )
00598          end_point1 = sub_loops_bounding_points1->next(0);
00599       else
00600          end_point1 = sub_loops_bounding_points1->next(jj+1);
00601       CubitVector uv1 = start_point1->coordinates();
00602       CubitVector uv2 = end_point1->coordinates();
00603 
00604         //Start with the current loop in the list.  The previous
00605         //lists have already been checked against this list.
00606       for ( kk = ii; kk < loops_bounding_points.size(); kk++ )
00607       {
00608         DLIList<CubitPoint*>* sub_loops_bounding_points2 = loops_bounding_points.next(kk);
00609         for ( ll = 0; ll < sub_loops_bounding_points2->size(); ll++ )
00610         {
00611           CubitPoint* start_point2 = sub_loops_bounding_points2->next(ll);
00612           CubitPoint* end_point2 = NULL;
00613           if ( ll == sub_loops_bounding_points2->size()-1 )
00614             end_point2 = sub_loops_bounding_points2->next(0);
00615           else
00616             end_point2 = sub_loops_bounding_points2->next(ll+1);
00617           CubitVector uv3 = start_point2->coordinates();
00618           CubitVector uv4 = end_point2->coordinates();
00619 
00620           if ( start_point1 == end_point2 || start_point1== start_point2 ||
00621                end_point1 == start_point2 || end_point1 == end_point2 )
00622             continue;
00623           
00624           if ( (uv1 == uv3) && (uv2 == uv4) )
00625             continue;
00626           
00627           double denom = ( (uv4.y() - uv3.y()) * (uv2.x() - uv1.x() ) -
00628                            (uv4.x() - uv3.x()) * (uv2.y() - uv1.y() ) );
00629           
00630           double numer_a = ( (uv4.x() - uv3.x()) * (uv1.y() - uv3.y() ) -
00631                              (uv4.y() - uv3.y()) * (uv1.x() - uv3.x() ) );
00632           
00633           double numer_b = ( (uv2.x() - uv1.x()) * (uv1.y() - uv3.y() ) -
00634                              (uv2.y() - uv1.y()) * (uv1.x() - uv3.x() ) );
00635           
00636           if( denom < CUBIT_RESABS && denom > -CUBIT_RESABS )
00637           {
00638               //The lines are parallel.  Check the numerators to
00639               //see if they are coincident.
00640             if( numer_a < CUBIT_RESABS && numer_a > -CUBIT_RESABS &&
00641                 numer_b < CUBIT_RESABS && numer_b > -CUBIT_RESABS )
00642             {
00643                 //The lines are coincident.  Now see if the segments
00644                 //over lap.
00645               double ua;
00646               double ub;
00647               if( (uv2.x() - uv1.x()) < CUBIT_RESABS &&
00648                   (uv2.x() - uv1.x()) > -CUBIT_RESABS )
00649               {
00650         if( (uv2.y() - uv1.y() < CUBIT_RESABS &&
00651              uv2.y() - uv1.y() > -CUBIT_RESABS)){
00652             PRINT_ERROR("Can't mesh due to zero-length edge.\n");
00653             return CUBIT_FAILURE;
00654         }
00655                 ua = (uv3.y() - uv1.y() ) / ( uv2.y() - uv1.y() );
00656                 ub = (uv4.y() - uv1.y() ) / ( uv2.y() - uv1.y() );
00657               }
00658               else
00659               {
00660                 ua = (uv3.x() - uv1.x() ) / ( uv2.x() - uv1.x() );
00661                 ub = (uv4.x() - uv1.x() ) / ( uv2.x() - uv1.x() );
00662               }
00663                   
00664               if ((ua > CUBIT_RESABS && ua < (1.0 - CUBIT_RESABS)) ||
00665                   (ub > CUBIT_RESABS && ub < (1.0 - CUBIT_RESABS)))
00666               {
00667                   //The segments overlap.
00668                 PRINT_ERROR("Can't mesh due to coincident edges \n\n\n");
00669                 return CUBIT_FAILURE;
00670               }
00671             }
00672           }
00673           else
00674           {
00675             double ua = numer_a / denom;
00676             double ub = numer_b / denom;
00677 
00678             if( ( (ua > CUBIT_RESABS && ua < (1.0-CUBIT_RESABS)) &&
00679                   (ub > CUBIT_RESABS && ub < (1.0-CUBIT_RESABS))) )
00680             {
00681               PRINT_ERROR("Can't mesh due to self intersecting edges \n\n\n");
00682               return CUBIT_FAILURE;
00683             }
00684           }
00685         }//end of for ( ll )
00686       }// end of for ( kk )
00687     }//end of for ( jj )
00688   }// end of for ( ii )
00689 
00690   return CUBIT_SUCCESS;
00691 }
00692 
00693 bool LoopParamTool::double_equal(double val, double equal_to)
00694 {
00695   if ( (val < (equal_to + CUBIT_RESABS)) && (val > (equal_to - CUBIT_RESABS)) )
00696     return true;
00697   else
00698     return false;
00699 }
00700 
00701 //EOF
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines