cgma
FBPolyhedron.cpp
Go to the documentation of this file.
00001 /*
00002  *
00003  *
00004  * Copyright (C) 2004 Sandia Corporation.  Under the terms of Contract DE-AC04-94AL85000
00005  * with Sandia Corporation, the U.S. Government retains certain rights in this software.
00006  *
00007  * This file is part of facetbool--contact via cubit@sandia.gov
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  * This library is distributed in the hope that it will be useful,
00015  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017  * Lesser General Public License for more details.
00018  *
00019  * You should have received a copy of the GNU Lesser General Public
00020  * License along with this library; if not, write to the Free Software
00021  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022  *
00023  *
00024  *
00025  */
00026 
00027 #include <math.h>
00028 #include <deque>
00029 #include <algorithm>
00030 
00031 #include "FBPolyhedron.hpp"
00032 #include "FBStructs.hpp"
00033 #include "FBRetriangulate.hpp"
00034 #include "CubitDefines.h"
00035 #include "CubitMessage.hpp"
00036 #include "CubitVector.hpp"
00037 #include "GfxDebug.hpp"
00038 const double BOX_CRACK = 1.e-4;
00039 
00040 FBPolyhedron::FBPolyhedron()
00041 {
00042 
00043   polyxmin = polyymin = polyzmin = CUBIT_DBL_MAX;
00044   polyxmax = polyymax = polyzmax = -polyxmin;
00045   hashobj = new IntegerHash(NUMHASHBINS,20);
00046   original_numtris = 0;
00047   
00048 }
00049 
00050 FBPolyhedron::~FBPolyhedron()
00051 {
00052   unsigned int i;
00053 
00054   delete hashobj; 
00055   delete kdtree;
00056   for ( i = 0; i < verts.size(); i++ ) {
00057     delete verts[i];
00058   } 
00059 
00060   for ( i = 0; i < tris.size(); i++ ) {
00061     delete tris[i];  
00062   }
00063 
00064   for ( i = 0; i < dudded_tris.size(); i++ ) {
00065     delete dudded_tris[i];
00066   }
00067 
00068   for ( i = 0; i < orphaned_edges.size(); i++ )
00069   {
00070     delete orphaned_edges[i];
00071   }
00072 }
00073 
00074 
00075 
00076 
00077 CubitStatus FBPolyhedron::makepoly(const std::vector<double>& coords,
00078                                    const std::vector<int>& connections,
00079                                    std::vector<int> *f_c_indices)
00080 {
00081   int hashvalue, parent, cubitfacetindex;
00082   int cubitedge0index, cubitedge1index, cubitedge2index;
00083   unsigned int i;
00084   FB_Coord *mycoord; 
00085   FB_Triangle *mytri;
00086   CubitStatus status;
00087   FSBOXVECTOR boxvector;
00088   std::vector<int>::iterator dpi;
00089   status = CUBIT_SUCCESS;
00090   
00091    for ( i = 0; i < coords.size(); i += 3 ) {
00092 
00093       mycoord = new FB_Coord(coords[i],coords[i+1],coords[i+2]); 
00094       hashvalue = makeahashvaluefrom_coord(coords[i],coords[i+1],coords[i+2]);
00095       hashobj->addtoHashList(hashvalue,verts.size());
00096       verts.push_back(mycoord);  
00097    }
00098    numpts = verts.size();
00099    parent = -1;
00100    if ( f_c_indices ){
00101      dpi = f_c_indices->begin();
00102    }
00103    for ( i = 0; i < connections.size(); i += 3 ) {
00104        //as a safety check make, not only do we see if we have the
00105        // f_c_indices vector, but make sure we don't go past the end of
00106        // that vector.
00107      if ( f_c_indices  && (i<f_c_indices->size())){
00108        cubitfacetindex = *dpi;
00109        cubitedge0index = *(dpi+1);
00110        cubitedge1index = *(dpi+2);
00111        cubitedge2index = *(dpi+3);
00112        dpi += 4;    
00113      }
00114      else {
00115        cubitfacetindex = 0;
00116        cubitedge0index = cubitedge1index = cubitedge2index = 0;
00117      }
00118      mytri = new FB_Triangle(connections[i],connections[i+1],connections[i+2],
00119                              parent,cubitfacetindex,cubitedge0index,
00120                              cubitedge1index,cubitedge2index);
00121      make_tri_boundingbox(mytri);
00122 
00123      boxvector.push_back(&mytri->boundingbox);  // for the KdTree
00124 
00125      make_tri_plane_coeffs(mytri);     
00126     
00127      tris.push_back(mytri);
00128      
00129    }
00130  
00131    numtris = original_numtris = tris.size();
00132    kdtree = new KDTree();
00133    kdtree->makeKDTree(numtris,boxvector);
00134    
00135    return status;        
00136 }
00137       
00138 bool FBPolyhedron::boxintersection(FBPolyhedron *otherpoly)
00139 {
00140 double xmin, ymin, zmin, xmax, ymax, zmax;
00141 
00142   xmin = otherpoly->polyxmin;
00143   ymin = otherpoly->polyymin;
00144   zmin = otherpoly->polyzmin;
00145   xmax = otherpoly->polyxmax;
00146   ymax = otherpoly->polyymax;
00147   zmax = otherpoly->polyzmax;
00148   
00149   if ( (polyxmin > xmax) || (polyxmax < xmin) ) return false;
00150   if ( (polyymin > ymax) || (polyymax < ymin) ) return false;
00151   if ( (polyzmin > zmax) || (polyzmax < zmin) ) return false;
00152    
00153   return true;
00154 }
00155 
00156 int FBPolyhedron::makeahashvaluefrom_coord(double x, double y, double z)
00157 {
00158 double mantissasum;
00159 
00160       if ( fabs(x) < 1.e-3 ) x = 0.0;
00161       if ( fabs(y) < 1.e-3 ) y = 0.0;
00162       if ( fabs(z) < 1.e-3 ) z = 0.0;
00163       mantissasum = (int)(10000.0*fabs(x) + 0.5) + 
00164                     (int)(10000.0*fabs(y) + 0.5) + 
00165           (int)(10000.0*fabs(z) + 0.5);
00166       
00167       return (int)(mantissasum) % NUMHASHBINS;
00168 }
00169 
00170 int FBPolyhedron::addavertex(double x, double y, double z)
00171 {
00172 int hashvalue, i, ifoundit;
00173 int *hasharrayptr, hasharraysize;
00174 double xval, yval, zval;
00175 FB_Coord *mycoord, *newcoord;
00176 
00177   hashvalue = makeahashvaluefrom_coord(x,y,z);
00178   
00179   hasharrayptr = hashobj->getHashBin(hashvalue,&hasharraysize);
00180   
00181   ifoundit = -1;
00182   for ( i = 0; i < hasharraysize; i++ ) {
00183     mycoord = verts[hasharrayptr[i]];
00184     xval = mycoord->coord[0];
00185     yval = mycoord->coord[1];
00186     zval = mycoord->coord[2];
00187     if ( ( fabs(xval-x) < EPSILON ) && 
00188     ( fabs(yval-y) < EPSILON ) &&
00189     ( fabs(zval-z) < EPSILON ) ) {
00190       ifoundit = hasharrayptr[i];
00191       break;
00192     }
00193   }
00194   if ( ifoundit == -1 ) {
00195     newcoord = new FB_Coord(x,y,z);
00196     ifoundit = verts.size();    
00197     verts.push_back(newcoord);
00198     hashobj->addtoHashList(hashvalue,ifoundit);
00199   }
00200   verts[ifoundit]->is_on_boundary = true;
00201   return ifoundit;
00202   
00203 }
00204 
00205 void FBPolyhedron::putnewpoints(std::vector<double>& newpoints)
00206 {
00207 double coordinate;
00208 unsigned int i;
00209 
00210 //  numpts is the original number of points.  Any new points were added  
00211 //  onto the end of the verts vector.
00212 
00213   if ( verts.size() > numpts ) {
00214     for ( i = numpts; i < verts.size(); i++ ) {
00215       coordinate = verts[i]->coord[0];
00216       newpoints.push_back(coordinate);
00217       coordinate = verts[i]->coord[1];
00218       newpoints.push_back(coordinate);
00219       coordinate = verts[i]->coord[2];
00220       newpoints.push_back(coordinate);
00221     }  
00222   }
00223 }
00224 
00225 void FBPolyhedron::putedges(std::vector<int>& newedges)
00226 {
00227 std::deque<unsigned int> edgedeque;
00228 std::multimap<unsigned int,unsigned int>::iterator p3;
00229 unsigned int startedge, startvert, edgnum;
00230 std::multimap<unsigned int,unsigned int>::iterator pub;
00231 unsigned int first, second;
00232 unsigned int i, numedgesadded;
00233 int v0, v1;
00234 bool foundone;
00235 
00236   for ( i = 0; i < intersection_edges.size(); i++ ) {
00237     v0 = intersection_edges[i]->v0;
00238     v1 = intersection_edges[i]->v1;
00239 
00240     edgnum = i;
00241     edgmmap.insert(std::pair<const unsigned int, unsigned int>(v0,edgnum));        
00242     edgmmap.insert(std::pair<const unsigned int, unsigned int>(v1,edgnum));    
00243   }
00244     numedgesadded = 0;
00245     startedge = 0;
00246 
00247   while ( numedgesadded < intersection_edges.size() ) {
00248     if ( intersection_edges[startedge]->mark == true ) {
00249       startedge++;
00250       continue;
00251     }
00252     startvert = v0 = intersection_edges[startedge]->v0;
00253     v1 = intersection_edges[startedge]->v1;
00254     intersection_edges[startedge]->mark = true;
00255     numedgesadded += 1;
00256     edgedeque.push_back(v0);
00257     edgedeque.push_back(v1);
00258  
00259     foundone = true;
00260 
00261     while ( foundone == true ) {       
00262       p3 = edgmmap.find(v1);
00263       pub = edgmmap.upper_bound(v1);
00264 
00265       if ( p3 != edgmmap.end() ) {
00266         do {
00267           foundone = false;
00268      first = p3->first;
00269      second = p3->second;
00270      if ( (second == startedge) && (first == startvert) ) {
00271        break;
00272      }
00273      if ( (intersection_edges[second]->v0 != v0) &&
00274           (intersection_edges[second]->v1 != v0) ) {
00275        if ( intersection_edges[second]->mark == true ) {
00276          p3++;
00277          continue;
00278        }
00279        if (intersection_edges[second]->v0 == v1 ) {
00280          v1 = intersection_edges[second]->v1;
00281          v0 = intersection_edges[second]->v0;
00282        } else {
00283          v1 = intersection_edges[second]->v0;
00284          v0 = intersection_edges[second]->v1;
00285        }
00286             intersection_edges[second]->mark = true;
00287             numedgesadded += 1;
00288             edgedeque.push_back(v1);
00289             foundone = true;
00290      }
00291      p3++;
00292         } while ( (p3 != pub) && (foundone == false) );
00293         if ( foundone == false ) break;
00294       } 
00295     } //                     end while ( foundone == true )
00296     v1 = startvert;
00297     foundone = true;
00298     while ( foundone == true ) {       
00299         p3 = edgmmap.find(v1);
00300         pub = edgmmap.upper_bound(v1);
00301 
00302         if ( p3 != edgmmap.end() ) {
00303           do {
00304             foundone = false;
00305        first = p3->first;
00306        second = p3->second;
00307        if ( (intersection_edges[second]->v0 == v1) ||
00308             (intersection_edges[second]->v1 == v1) ) {
00309          if ( intersection_edges[second]->mark == true ) {
00310            p3++;
00311            continue;
00312          }
00313          if (intersection_edges[second]->v0 == v1 ) {
00314            v1 = intersection_edges[second]->v1;
00315            v0 = intersection_edges[second]->v0;
00316          } else {
00317            v1 = intersection_edges[second]->v0;
00318            v0 = intersection_edges[second]->v1;
00319          }
00320               intersection_edges[second]->mark = true;
00321          numedgesadded += 1;
00322               edgedeque.push_front(v1);
00323          foundone = true;
00324        }
00325        p3++;
00326           } while ( (p3 != pub) && (foundone == false) );
00327           if ( foundone == false ) break;
00328         } 
00329     }
00330   
00331     unsigned int size = edgedeque.size();
00332 
00333     if ( size > 0 ) {
00334       newedges.push_back(size);
00335       for ( i = 0; i < size; i++ ) {
00336         newedges.push_back(edgedeque[i]);
00337       }
00338     }
00339     edgedeque.clear();
00340     startedge++;
00341   } //                     end while ( numedgesadded < intersection_edges.size() )
00342  
00343   newedges.push_back(0);
00344 
00345 
00346 }
00347 
00348 bool FBPolyhedron::edge_exists(int v0, int v1) 
00349 {
00350 bool exists = false;
00351 unsigned int i;
00352 
00353   for ( i = 0; i < intersection_edges.size(); i++ ) {
00354     if ( ( (intersection_edges[i]->v0 == v0) && 
00355            (intersection_edges[i]->v1 == v1) ) ||
00356          ( (intersection_edges[i]->v0 == v1) && 
00357            (intersection_edges[i]->v1 == v0) ) ) {
00358       exists = true;
00359       break;
00360     }   
00361   }
00362 
00363   return exists;
00364 }
00365 
00366 CubitStatus FBPolyhedron::retriangulate(std::vector<int>& newfacets, 
00367                                         std::vector<int>& newfacetsindex)
00368 {
00369 CubitStatus status;
00370 FBRetriangulate *retriangulater;
00371 unsigned int i;
00372 
00373   status = CUBIT_SUCCESS;
00374   for ( i = 0; i < tris.size(); i++ ) {
00375     if ( tris[i]->dudded == true ) {
00376       tris[i]->parent = (int)i;
00377       retriangulater = new FBRetriangulate(verts, tris, newfacets, newfacetsindex); 
00378       status = retriangulater->retriangulate_this_tri(i, orphaned_edges);
00379       delete retriangulater;
00380     }
00381   
00382   }
00383   
00384   return status;
00385 }
00386 
00387 CubitStatus FBPolyhedron::retriangulate(std::vector<int>& newfacets)
00388 {
00389 CubitStatus status;
00390 FBRetriangulate *retriangulater;
00391 unsigned int i;
00392 
00393   status = CUBIT_SUCCESS;
00394   for ( i = 0; i < tris.size(); i++ ) {
00395     if ( tris[i]->dudded == true ) {
00396       tris[i]->parent = (int)i;
00397       retriangulater = new FBRetriangulate(verts, tris, newfacets); 
00398       status = retriangulater->retriangulate_this_tri(i, orphaned_edges);
00399       std::vector<FB_Edge*>::iterator iter;
00400       delete retriangulater;
00401     }
00402   
00403   }
00404   
00405   return status;
00406 }
00407 
00408 bool FBPolyhedron::edge_exists_in_tri(FB_Triangle& tri, int v0, int v1)
00409 {
00410 FB_Edge *edge;
00411 std::vector<FB_Edge*>::iterator dp;
00412 
00413   dp = tri.edge_list.begin();
00414   while ( dp != tri.edge_list.end() ) {
00415     edge = *dp;
00416     if ( ( ( (edge->v0) == v0 ) && ( (edge->v1) == v1 ) ) ||
00417          ( ( (edge->v0) == v1 ) && ( (edge->v1) == v0 ) ) )
00418       return true;
00419     dp++;  
00420   }
00421   
00422   return false;
00423 }
00424 
00425 void FBPolyhedron::add_new_triangle_data()
00426 {
00427 unsigned int i;
00428 FB_Triangle *tri;
00429 
00430   for ( i = original_numtris; i < tris.size(); i++ ) {
00431   
00432     tri = tris[i];
00433 //  make the bounding box
00434     make_tri_boundingbox(tri);
00435 //  make the plane coefficients
00436     make_tri_plane_coeffs(tri);
00437   }
00438   
00439 }
00440 
00441 void FBPolyhedron::make_tri_plane_coeffs(FB_Triangle *tri)
00442 {
00443 FB_Coord *mycoord;
00444 double x1, x2, x3, y1, y2, y3, z1, z2, z3, e1x, e1y, e1z, e2x, e2y, e2z;
00445 double a, b, c, d, dtemp;
00446 
00447      mycoord = verts[tri->v0];
00448      x1 = mycoord->coord[0];
00449      y1 = mycoord->coord[1];
00450      z1 = mycoord->coord[2];
00451      mycoord = verts[tri->v1];
00452      x2 = mycoord->coord[0];
00453      y2 = mycoord->coord[1];
00454      z2 = mycoord->coord[2];
00455      mycoord = verts[tri->v2];
00456      x3 = mycoord->coord[0];
00457      y3 = mycoord->coord[1];
00458      z3 = mycoord->coord[2];
00459      e1x = x1 - x2; e1y = y1 - y2; e1z = z1 - z2;
00460      e2x = x3 - x2; e2y = y3 - y2; e2z = z3 - z2;
00461      a = e1z*e2y - e2z*e1y;
00462      b = e1x*e2z - e2x*e1z;
00463      c = e1y*e2x - e2y*e1x;
00464      dtemp = sqrt(a*a + b*b + c*c);
00465      if ( dtemp > EPSILON2 ) {
00466        a /= dtemp;
00467        b /= dtemp;
00468        c /= dtemp;
00469      } else {
00470        PRINT_WARNING("small-area triangle\n");
00471      }
00472      d = -(a*x1 + b*y1 + c*z1);
00473      tri->a = a; tri->b = b; tri->c = c; tri->d = d;
00474 
00475 }
00476 
00477 void FBPolyhedron::make_tri_boundingbox(FB_Triangle *tri)
00478 {
00479 double xmin, ymin, zmin, xmax, ymax, zmax;
00480 int j;
00481 int connections[3];
00482 FB_Coord *mycoord;
00483 
00484      xmin = ymin = zmin = CUBIT_DBL_MAX;
00485      xmax = ymax = zmax = -xmin;
00486      connections[0] = tri->v0; connections[1] = tri->v1; connections[2] = tri->v2;
00487      for ( j = 0; j < 3; j++ ) { // make the bounding box
00488        mycoord = verts[connections[j]];
00489        xmin = ( xmin < mycoord->coord[0] ) ? xmin : mycoord->coord[0];
00490        xmax = ( xmax > mycoord->coord[0] ) ? xmax : mycoord->coord[0];
00491        ymin = ( ymin < mycoord->coord[1] ) ? ymin : mycoord->coord[1];
00492        ymax = ( ymax > mycoord->coord[1] ) ? ymax : mycoord->coord[1];
00493        zmin = ( zmin < mycoord->coord[2] ) ? zmin : mycoord->coord[2];
00494        zmax = ( zmax > mycoord->coord[2] ) ? zmax : mycoord->coord[2];
00495      }
00496      if ( (xmax - xmin) < BOX_CRACK ) {
00497        xmax += BOX_CRACK; xmin -= BOX_CRACK;
00498      }
00499      if ( (ymax - ymin) < BOX_CRACK ) {
00500        ymax += BOX_CRACK; ymin -= BOX_CRACK;
00501      }
00502      if ( (zmax - zmin) < BOX_CRACK ) {
00503        zmax += BOX_CRACK; zmin -= BOX_CRACK;
00504      }     
00505      tri->boundingbox.xmin = xmin; tri->boundingbox.xmax = xmax;
00506      tri->boundingbox.ymin = ymin; tri->boundingbox.ymax = ymax;
00507      tri->boundingbox.zmin = zmin; tri->boundingbox.zmax = zmax;
00508 
00509      //  Update the object's bounding box
00510      polyxmin = ( polyxmin < xmin ) ? polyxmin : xmin;
00511      polyymin = ( polyymin < ymin ) ? polyymin : ymin;
00512      polyzmin = ( polyzmin < zmin ) ? polyzmin : zmin;
00513      polyxmax = ( polyxmax > xmax ) ? polyxmax : xmax;
00514      polyymax = ( polyymax > ymax ) ? polyymax : ymax;
00515      polyzmax = ( polyzmax > zmax ) ? polyzmax : zmax;
00516 
00517 }
00518 
00519 void FBPolyhedron::removeduddedtriangles()  //  Compacts the tris vector
00520 {
00521 unsigned int i;
00522 int j;
00523 
00524   j = 0;
00525   for ( i = 0; i < tris.size(); i++ ) {
00526     if ( tris[i]->dudded == true )
00527     {
00528       dudded_tris.push_back(tris[i]);
00529     }
00530     else
00531     {
00532       tris[j] = tris[i];
00533       j++;
00534     }
00535   }
00536   tris.resize(j);
00537 }
00538 
00539   //find the largest and smallest angles in this triangle
00540 bool FBPolyhedron::min_max_angles_in_fb_triangle(FB_Triangle *triangle,
00541                                                 double& min_angle,
00542                                                 double& max_angle)
00543 {
00544   CubitVector vert_0(verts[triangle->v0]->coord[0],
00545                      verts[triangle->v0]->coord[1],
00546                      verts[triangle->v0]->coord[2]);
00547   CubitVector vert_1(verts[triangle->v1]->coord[0],
00548                      verts[triangle->v1]->coord[1],
00549                      verts[triangle->v1]->coord[2]);
00550   CubitVector vert_2(verts[triangle->v2]->coord[0],
00551                      verts[triangle->v2]->coord[1],
00552                      verts[triangle->v2]->coord[2]);
00553   CubitVector sides[3];
00554   sides[0] = vert_1 - vert_0;
00555   sides[1] = vert_2 - vert_1;
00556   sides[2]= vert_0 - vert_2;
00557   if(sides[0].length_squared() < EPSILON ||
00558      sides[1].length_squared() < EPSILON ||
00559      sides[2].length_squared() < EPSILON ){
00560     min_angle =0.0;
00561     max_angle =180.0;
00562     return false;
00563   }
00564   double curr_angle;
00565   min_angle = sides[1].interior_angle(-sides[0]);
00566   max_angle = min_angle;
00567   curr_angle = sides[2].interior_angle(-sides[1]);
00568   if(curr_angle<min_angle){
00569     min_angle=curr_angle;
00570   }
00571   if(curr_angle>max_angle){
00572     min_angle=curr_angle;
00573   }
00574   curr_angle = sides[0].interior_angle(-sides[2]);
00575   if(curr_angle<min_angle){
00576     min_angle=curr_angle;
00577   }
00578   if(curr_angle>max_angle){
00579     min_angle=curr_angle;
00580   }
00581   return true;
00582   
00583 }
00584 //draw the "boundary edges" of the polyhedron.
00585 void FBPolyhedron::debug_draw_boundary_edges(int color)
00586 {
00587   unsigned int i;
00588   FB_Triangle* triangle=NULL;
00589   for(i = 0; i<tris.size(); i++){
00590     if(!tris[i]->dudded){
00591       triangle = tris[i];
00592       unsigned int counter = 0;
00593       CubitVector vert_0(verts[triangle->v0]->coord[0],
00594                            verts[triangle->v0]->coord[1],
00595                          verts[triangle->v0]->coord[2]);
00596       CubitVector vert_1(verts[triangle->v1]->coord[0],
00597                          verts[triangle->v1]->coord[1],
00598                          verts[triangle->v1]->coord[2]);
00599       CubitVector vert_2(verts[triangle->v2]->coord[0],
00600                          verts[triangle->v2]->coord[1],
00601                          verts[triangle->v2]->coord[2]);
00602       if(triangle->cubitedge0index){
00603         counter++;
00604         GfxDebug::draw_line(vert_0, vert_1, color);
00605       }
00606       if(triangle->cubitedge1index){
00607         counter++;
00608         GfxDebug::draw_line(vert_1, vert_2, color);
00609       }
00610       if(triangle->cubitedge2index){
00611         counter++;
00612         GfxDebug::draw_line(vert_2, vert_0, color);
00613       }
00614       if(counter != triangle->edge_list.size()){
00615         PRINT_WARNING("Possible debug problem.\n");
00616       }
00617     }
00618   }
00619 }
00620 
00621 //draw a single triangle
00622 void FBPolyhedron::debug_draw_fb_triangle(FB_Triangle *triangle)
00623 {
00624   CubitVector vert_0(verts[triangle->v0]->coord[0],
00625                      verts[triangle->v0]->coord[1],
00626                      verts[triangle->v0]->coord[2]);
00627   CubitVector vert_1(verts[triangle->v1]->coord[0],
00628                      verts[triangle->v1]->coord[1],
00629                      verts[triangle->v1]->coord[2]);
00630   CubitVector vert_2(verts[triangle->v2]->coord[0],
00631                      verts[triangle->v2]->coord[1],
00632                      verts[triangle->v2]->coord[2]);
00633   GfxDebug::draw_point(vert_0, CUBIT_RED_INDEX);
00634   GfxDebug::draw_point(vert_1, CUBIT_RED_INDEX);
00635   GfxDebug::draw_point(vert_2, CUBIT_RED_INDEX);
00636   GfxDebug::draw_line(vert_0, vert_1, CUBIT_BLUE_INDEX);
00637   GfxDebug::draw_line(vert_1, vert_2, CUBIT_BLUE_INDEX);
00638   GfxDebug::draw_line(vert_2, vert_0, CUBIT_BLUE_INDEX);
00639   int draw_normal = 1;
00640   if(draw_normal){
00641     CubitVector center_pos = (vert_0+vert_1+vert_2)/3.0;
00642     CubitVector opp_pos = center_pos + (vert_1-vert_0)*(vert_2-vert_0);
00643     GfxDebug::draw_line(center_pos, opp_pos, CUBIT_WHITE_INDEX);
00644   }
00645 }
00646 
00647 //get the largest and smallest angles in the polyhedron.  That is,
00648 //find the largest angle in any triangle in the polyhedron and find
00649 //the smallest angle in any triangle in the polyhedron.
00650 bool FBPolyhedron::min_max_angles_in_polyhedron(double& min_angle,
00651                                                   double& max_angle)
00652 {
00653   double curr_min;
00654   double curr_max;
00655   if(tris.size() < 1){
00656     min_angle = 0.0;
00657     max_angle = 180.0;
00658     return false;
00659   }
00660   min_max_angles_in_fb_triangle(tris[0], min_angle, max_angle);
00661   unsigned int i;
00662   for(i = 1; i<tris.size(); i++){
00663     if(!min_max_angles_in_fb_triangle(tris[i], curr_min, curr_max)){
00664       return false;
00665     }
00666     if(curr_min<min_angle){
00667       min_angle=curr_min;
00668     }
00669     if(curr_max>max_angle){
00670       max_angle=curr_max;
00671     }
00672   }
00673   return true;
00674 }
00675 
00676 #ifdef KEEP_BOYD10_KEEP
00677 //function that attempts to remove small angles be flipping edges...
00678 bool FBPolyhedron::remove_small_angles_in_triangle_range(int lower_index,
00679                                                          int upper_index)
00680 {
00681     //we currently only handle obtuse, degenerate triangles...
00682     //see if the triangle falls in this range
00683   int mydebug = 0;
00684   const double min_angle_in_degrees = 1;
00685   const double max_angle_in_degrees = 120;
00686   int j;
00687   for(j=lower_index;j<upper_index;j++){
00688       //this triangle has an extremely small angl
00689     double min_angle, max_angle;
00690     if(!min_max_angles_in_fb_triangle(tris[j], min_angle, max_angle))
00691       return false;
00692     
00693     if(min_angle < min_angle_in_degrees && max_angle>max_angle_in_degrees)
00694     {
00695         //PRINT_INFO("Found a degenerate, obtuse triangle.\n");
00696         //two vertices on longest edge of degenerate triangle
00697       int longest_edge_v1 = 0;
00698       int longest_edge_v2 = 0;
00699         //all vertices on triangle
00700       int v_indices[3];
00701       v_indices[0]=tris[j]->v0;
00702       v_indices[1]=tris[j]->v1;
00703       v_indices[2]=tris[j]->v2;
00704       double longest_length = -1.0;
00705       double temp_length = 0.0;
00706       int k=0;
00707         //find the longest edge
00708       for(k=0;k<3;k++){
00709         CubitVector v1(verts[v_indices[k]]->coord[0],
00710                        verts[v_indices[k]]->coord[1],
00711                        verts[v_indices[k]]->coord[2]);
00712         CubitVector v2(verts[v_indices[(k+1)%3]]->coord[0],
00713                        verts[v_indices[(k+1)%3]]->coord[1],
00714                        verts[v_indices[(k+1)%3]]->coord[2]);
00715         temp_length = (v1-v2).length_squared();
00716               
00717         if(temp_length>longest_length){
00718           longest_length = temp_length;
00719           longest_edge_v1 = v_indices[k];
00720           longest_edge_v2 = v_indices[(k+1)%3];
00721         }
00722       }
00723         //look for the other triangle in this range that shares this edge
00724       int found_it = 0;
00725       k=0;
00726       while (!found_it && k < (int)tris.size()){
00727         if(k!=j && tris[k]->are_fb_coords_in_fb_triangle(
00728           longest_edge_v1,longest_edge_v2)){
00729           found_it=1;
00730         }
00731         if(!found_it){
00732           k++;
00733         }
00734       }
00735         //we have found the two triangles adjacent to the longest
00736         //edge of the degenerate tri.
00737       if(!found_it){
00738           //PRINT_WARNING("Didn't find the other triangle\n");
00739         return false;
00740       }
00741 //       PRINT_INFO("Found it... j = %i, k = %i\n",j,k);
00742 //       PRINT_INFO("edge lists sizes %i, %i\n",tris[j]->edge_list.size(),tris[k]->edge_list.size());
00743       double min_angle_before_swap;
00744       min_max_angles_in_fb_triangle(tris[j], min_angle_before_swap,
00745                                     max_angle);
00746       // PRINT_INFO("Inititial tri j angle = %f\n",min_angle);
00747       min_max_angles_in_fb_triangle(tris[k], min_angle, max_angle);
00748       if(min_angle<min_angle_before_swap){
00749         min_angle_before_swap = min_angle;
00750       }
00751       if(mydebug){
00752         PRINT_INFO("Inititial tri k angle = %f\n",min_angle);
00753         GfxDebug::clear();
00754         debug_draw_fb_triangle(tris[j]);
00755         debug_draw_fb_triangle(tris[k]);    
00756         GfxDebug::mouse_xforms();
00757       }
00758         //now do a swap to remove the longest edge...
00759       int other_vj;
00760       int other_vk;
00761       int index_for_j = 0;
00762       int index_for_k =0;
00763       int old_value_for_j=0;
00764       int old_value_for_k=0;
00765       if(tris[j]->v0 == longest_edge_v1 ||
00766          tris[j]->v0 == longest_edge_v2 ){
00767         if(tris[j]->v1 == longest_edge_v1 ||
00768            tris[j]->v1 == longest_edge_v2){
00769           other_vj =tris[j]->v2;
00770           index_for_j=1;
00771         }
00772         else{
00773           other_vj =tris[j]->v1;
00774           index_for_j=0;
00775         }
00776       }
00777       else{
00778         other_vj =tris[j]->v0;
00779         index_for_j=2;
00780       }
00781       if(tris[k]->v0 == longest_edge_v1 ||
00782          tris[k]->v0 == longest_edge_v2 ){
00783         if(tris[k]->v1 == longest_edge_v1 ||
00784            tris[k]->v1 == longest_edge_v2){
00785           other_vk = tris[k]->v2;
00786           old_value_for_k=tris[k]->v1;
00787           tris[k]->v1=other_vj;
00788           index_for_k = 1;
00789         }
00790         else{
00791           other_vk = tris[k]->v1;
00792           old_value_for_k=tris[k]->v0;
00793           tris[k]->v0=other_vj;
00794           index_for_k=0;
00795         }
00796       }
00797       else{
00798         other_vk = tris[k]->v0;
00799         old_value_for_k=tris[k]->v2;
00800         tris[k]->v2=other_vj;
00801         index_for_k=2;
00802       }
00803       if(index_for_j == 0){
00804         old_value_for_j=tris[j]->v0;
00805         tris[j]->v0=other_vk;
00806       }
00807       else if(index_for_j == 1){
00808         old_value_for_j=tris[j]->v1;
00809         tris[j]->v1=other_vk;
00810       }
00811       else{
00812         old_value_for_j=tris[j]->v2;
00813         tris[j]->v2=other_vk;
00814       }
00815       double min_angle_after_swap;
00816       min_max_angles_in_fb_triangle(tris[j], min_angle_after_swap,
00817                                     max_angle);
00818         // PRINT_INFO("Modified tri j angle = %f\n",min_angle);
00819       min_max_angles_in_fb_triangle(tris[k], min_angle, max_angle);
00820       if(min_angle<min_angle_after_swap){
00821         min_angle_after_swap=min_angle;
00822       }
00823       if(min_angle_after_swap < min_angle_before_swap){
00824           // PRINT_WARNING("Swap was unsuccessful in removing bad angle.");
00825         if(index_for_j == 0){
00826           tris[j]->v0=old_value_for_j;
00827         }
00828         else if(index_for_j == 1){
00829           
00830           tris[j]->v1=old_value_for_j;
00831         }
00832         else{
00833           tris[j]->v2=old_value_for_j;
00834         }
00835         if(index_for_k == 0){
00836           tris[k]->v0=old_value_for_k;
00837         }
00838         else if(index_for_k == 1){
00839           
00840           tris[k]->v1=old_value_for_k;
00841         }
00842         else{
00843           tris[k]->v2=old_value_for_k;
00844         }
00845       }
00846       if(mydebug){
00847         PRINT_INFO("Modified tri k angle = %f\n",min_angle);
00848         GfxDebug::clear();
00849         debug_draw_fb_triangle(tris[j]);
00850         debug_draw_fb_triangle(tris[k]);    
00851         GfxDebug::mouse_xforms();
00852       }
00853     }
00854   }
00855   return true;
00856   
00857 }
00858 #endif
00859 
00860   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines