cgma
FBRetriangulate.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 <algorithm>
00029 #include "FBRetriangulate.hpp"
00030 #include "FBTiler.hpp"
00031 #include "CubitMessage.hpp"
00032 
00033 #ifdef KEEP_BOYD10_KEEP
00034 bool edgecf_less(const FB_Edge* efirst, const FB_Edge* esecond)
00035 {
00036   if ( efirst->quadrant < esecond->quadrant ) return true;
00037   if ( efirst->quadrant > esecond->quadrant ) return false;
00038 
00039   if ( efirst->slope < esecond->slope ) return true;
00040   else return false;       
00041 
00042 }
00043 #endif
00044     
00045 FBRetriangulate::FBRetriangulate(std::vector<FB_Coord *>& my_verts,
00046                                  std::vector<FB_Triangle *>& my_tris,
00047                                  std::vector<int>& my_newfacets,
00048                                  std::vector<int>& my_newfacetsindex)
00049 {
00050   verts = my_verts;
00051   tris = &my_tris;
00052   newfacets = &my_newfacets;
00053   newfacetsindex = &my_newfacetsindex;
00054   p_dir=0;
00055   s_dir=0;
00056 }
00057     
00058 FBRetriangulate::FBRetriangulate(std::vector<FB_Coord *>& my_verts,
00059                                  std::vector<FB_Triangle *>& my_tris,
00060                                  std::vector<int>& my_newfacets)
00061 {
00062   verts = my_verts;
00063   tris = &my_tris;
00064   newfacets = &my_newfacets;
00065   newfacetsindex = 0;
00066   p_dir=0;
00067   s_dir=0;
00068 }
00069 
00070 FBRetriangulate::~FBRetriangulate()
00071 {
00072 
00073 }
00074 
00075 CubitStatus FBRetriangulate::retriangulate_this_tri(int sequence, std::vector<FB_Edge*> &orphaned_edges)
00076 {
00077   CubitStatus status;
00078   double xspan, yspan, zspan;
00079   std::vector<FB_Triangle *>::iterator itt;
00080 
00081   itt = tris->begin();
00082   itt += sequence;
00083   my_tri = *itt;   
00084   status = CUBIT_SUCCESS;
00085   
00086 //  tri = &itt;  
00087   xspan = my_tri->boundingbox.xmax - my_tri->boundingbox.xmin;
00088   yspan = my_tri->boundingbox.ymax - my_tri->boundingbox.ymin;
00089   zspan = my_tri->boundingbox.zmax - my_tri->boundingbox.zmin;
00090   if ( ( fabs(my_tri->c) >= fabs(my_tri->a) ) &&
00091        ( fabs(my_tri->c) >= fabs(my_tri->b) ) ) {
00092     p_dir = 1; 
00093     s_dir = 0;
00094     if ( xspan > yspan ) {
00095       p_dir = 0; 
00096       s_dir = 1;
00097     }
00098   }
00099   else if ( fabs(my_tri->b) >= fabs(my_tri->a) )  {
00100     p_dir = 2; 
00101     s_dir = 0;
00102     if ( xspan > zspan ) {
00103       p_dir = 0;
00104       s_dir = 2;
00105     }
00106   }
00107   else {
00108     p_dir = 1;
00109     s_dir = 2;
00110     if ( zspan > yspan ) {
00111       p_dir = 2;
00112       s_dir = 1;
00113     }
00114   }
00115    
00116   if ( p_dir == 0 ) {
00117     if ( s_dir == 1 ) {
00118       if ( my_tri->c > 0.0 )
00119         winding = CCW;
00120       else
00121         winding = CW;
00122     }
00123     else if (s_dir == 2 ) {
00124       if ( my_tri->b > 0.0 )
00125         winding = CW;
00126       else
00127         winding = CCW;
00128     }
00129   }
00130   else if ( p_dir == 1 ) {
00131     if ( s_dir == 0 ) {
00132       if ( my_tri->c > 0.0 )
00133         winding = CW;
00134       else
00135         winding = CCW;
00136     }
00137     else if (s_dir == 2 ) {
00138       if ( my_tri->a > 0.0 )
00139         winding = CW;
00140       else
00141         winding = CCW;
00142     }
00143   }
00144   else if ( p_dir == 2 ) {
00145     if ( s_dir == 0 ) {
00146       if ( my_tri->b > 0.0 )
00147         winding = CCW;
00148       else
00149         winding = CW;
00150     }
00151     else if (s_dir == 1 ) {
00152       if ( my_tri->a > 0.0 )
00153         winding = CCW;
00154       else
00155         winding = CW;
00156     } 
00157   }
00158   else{
00159     PRINT_ERROR("Unexpected result.\n");
00160     return CUBIT_FAILURE;
00161   }
00162   
00163 //if ( winding == CCW ) winding = CW;
00164 //else winding = CCW;
00165   classify_edges();
00166   
00167 //  Add edges around the perimeter of the triangle to its edge list.
00168   status = add_bdry_edges(orphaned_edges);
00169   make_vert_list();
00170   
00171 //  Add internal edges where there is a local min or max 
00172   status = remove_min_max();
00173 
00174 // Get the vertex chains that will be retriangulated.  
00175   
00176   bool chain_status = true; 
00177   sort_vertstufflist_edges();  
00178   std::vector<int> *chainlist = new std::vector<int>;
00179 //  std::vector<FB_Triangle* > new_tris;
00180   unsigned int i, nfsize;
00181   nfsize = newfacets->size();
00182   while (chain_status == true ) {
00183     chain_status = get_a_chain(chainlist);
00184     if ( chain_status == false ) break;
00185     FBTiler *tiler = new FBTiler(verts,p_dir,s_dir,sequence,
00186                                  my_tri->a,my_tri->b,my_tri->c,newfacets);
00187     status = tiler->Tile_Poly(chainlist);
00188     delete tiler;
00189     chainlist->clear();
00190 
00191   }
00192 
00193   delete chainlist;
00194   unsigned int number_new;
00195   int e0index, e1index, e2index;
00196 
00197   number_new = newfacets->size() - nfsize;
00198   if ( number_new > 0 ) {
00199     std::vector<int>::iterator itp;
00200     itp = newfacets->begin();
00201     itp += nfsize;
00202     for ( i = 0; i < number_new; i += 3 ) {
00203       FB_Triangle *new_tri;
00204       e0index = e1index = e2index = 0;
00205       get_edge_indices(*itp,*(itp+1),*(itp+2),sequence,
00206                        e0index, e1index, e2index);
00207       new_tri = new FB_Triangle(*itp,*(itp+1),*(itp+2),
00208                                 sequence,my_tri->cubitsurfaceindex,
00209                                 e0index, e1index, e2index);
00210       tris->push_back(new_tri);                         
00211       itp += 3;   
00212     }
00213     if ( newfacetsindex ) newfacetsindex->push_back(nfsize);
00214   }  
00215   for ( i = 0; i < vertstufflist.size(); i++ ) {
00216     delete vertstufflist[i];
00217   }
00218 
00219   vertstufflist.clear();
00220 
00221   return status;
00222 
00223 }
00224 
00225 CubitStatus FBRetriangulate::remove_min_max()
00226 {
00227   CubitStatus status;
00228   FB_Edge *edge;
00229   std::vector<FB_Edge*>::iterator dpe;
00230   bool goes_down, goes_up;
00231   VertexStuff *vstuff;
00232   unsigned int i;
00233   int this_vert, that_vert;
00234   double this_vert_p_dir_coord, that_vert_p_dir_coord;
00235   
00236   status = CUBIT_SUCCESS;
00237 
00238   for ( i = 0; i < vertstufflist.size(); i++ ) {
00239     vstuff = vertstufflist[i];
00240     if ( vstuff->v0type != INTERIOR_VERT ) continue;
00241     goes_down = goes_up = false;
00242     this_vert = vstuff->v0;
00243     this_vert_p_dir_coord = verts[this_vert]->coord[p_dir];
00244     dpe = vstuff->edge_list.begin();
00245     while ( dpe != vstuff->edge_list.end() ) {
00246       edge = *dpe;
00247       dpe++;
00248       if ( edge->v0 == this_vert ) that_vert = edge->v1;
00249       else that_vert = edge->v0;
00250       that_vert_p_dir_coord = verts[that_vert]->coord[p_dir]; 
00251       if ( that_vert_p_dir_coord < this_vert_p_dir_coord )
00252         goes_down = true;
00253       else if ( that_vert_p_dir_coord > this_vert_p_dir_coord )
00254         goes_up = true;
00255       else {
00256         double this_vert_s_dir_coord, that_vert_s_dir_coord;
00257         this_vert_s_dir_coord = verts[this_vert]->coord[s_dir];
00258         that_vert_s_dir_coord = verts[that_vert]->coord[s_dir]; 
00259         if ( this_vert_s_dir_coord < that_vert_s_dir_coord )
00260           goes_up = true;
00261         else goes_down = true;
00262       }
00263     }
00264 
00265     if ( goes_down == false ) {
00266       add_edge_down(this_vert,i);
00267     } else if ( goes_up == false ) {
00268       add_edge_up(this_vert,i);    
00269     } 
00270   }
00271          
00272   return status;
00273   
00274 }
00275 
00276 CubitStatus FBRetriangulate::add_bdry_edges(std::vector<FB_Edge*> &orphaned_edges)
00277 {
00278 //  For each triangle edge, make a list of pairs of edge points and
00279 //  vertex numbers that touch the edge.  Then sort the list by distance
00280 //  from the start of the edge.  (This was what we really put in the
00281 //  pair, not the point itself.)  Finally, go through this sorted list
00282 //  and add edges to the edge_list.  If the edge_list is empty, we
00283 //  still must add in the triangle edge itself.
00284 
00285   CubitStatus status;
00286   FB_Edge *edge;
00287   std::list< std::pair<double,int> > edge0_list, edge1_list, edge2_list;
00288   std::list< std::pair<double,int> >::iterator dp;
00289   std::vector<FB_Edge*>::iterator dpe, dpe_orig;
00290   double vx0, vy0, vz0, vx1, vy1, vz1, vx2, vy2, vz2, dist;
00291   std::pair<double,int> mypair;
00292 
00293   status = CUBIT_SUCCESS;
00294 
00295   vx0 = verts[my_tri->v0]->coord[0];
00296   vy0 = verts[my_tri->v0]->coord[1];
00297   vz0 = verts[my_tri->v0]->coord[2];
00298   vx1 = verts[my_tri->v1]->coord[0];
00299   vy1 = verts[my_tri->v1]->coord[1];
00300   vz1 = verts[my_tri->v1]->coord[2];
00301   vx2 = verts[my_tri->v2]->coord[0];
00302   vy2 = verts[my_tri->v2]->coord[1];
00303   vz2 = verts[my_tri->v2]->coord[2];
00304 
00305   dpe = my_tri->edge_list.begin();
00306   while ( dpe != my_tri->edge_list.end() ) {
00307     edge = *dpe;
00308     dpe++;
00309     if ( (edge->v0_type == INTERIOR_VERT) &&
00310          (edge->v1_type == INTERIOR_VERT) )
00311       continue;
00312     if ( edge->v0_type == EDGE_0 ) {
00313       dist = get_dist(vx0,vy0,vz0,verts[edge->v0]->coord[0],
00314                       verts[edge->v0]->coord[1],verts[edge->v0]->coord[2]);
00315       edge0_list.push_back(std::pair<double,int>(dist,edge->v0));
00316     } else if ( edge->v0_type == EDGE_1 ) { 
00317       dist = get_dist(vx1,vy1,vz1,verts[edge->v0]->coord[0],
00318                       verts[edge->v0]->coord[1],verts[edge->v0]->coord[2]);
00319       edge1_list.push_back(std::pair<double,int>(dist,edge->v0));
00320     } else if ( edge->v0_type == EDGE_2 ) {
00321       dist = get_dist(vx2,vy2,vz2,verts[edge->v0]->coord[0],
00322                       verts[edge->v0]->coord[1],verts[edge->v0]->coord[2]);
00323       edge2_list.push_back(std::pair<double,int>(dist,edge->v0));
00324     } 
00325     if ( edge->v1_type == EDGE_0 ) {
00326       dist = get_dist(vx0,vy0,vz0,verts[edge->v1]->coord[0],
00327                       verts[edge->v1]->coord[1],verts[edge->v1]->coord[2]);
00328       edge0_list.push_back(std::pair<double,int>(dist,edge->v1));
00329     } else if ( edge->v1_type == EDGE_1 ) {
00330       dist = get_dist(vx1,vy1,vz1,verts[edge->v1]->coord[0],
00331                       verts[edge->v1]->coord[1],verts[edge->v1]->coord[2]);
00332       edge1_list.push_back(std::pair<double,int>(dist,edge->v1));
00333     } else if ( edge->v1_type == EDGE_2 ) {
00334       dist = get_dist(vx2,vy2,vz2,verts[edge->v1]->coord[0],
00335                       verts[edge->v1]->coord[1],verts[edge->v1]->coord[2]);
00336       edge2_list.push_back(std::pair<double,int>(dist,edge->v1));
00337     }     
00338   }
00339   
00340   edge0_list.sort(); 
00341   edge1_list.sort(); 
00342   edge2_list.sort();
00343    
00344 //  Now we have to remove all BDRY_EDGEs because they will be made anew
00345 //  in what follows.  Erasing elements from a vector is inefficient, but
00346 //  this shouldn't happen often.
00347   dpe = my_tri->edge_list.begin();
00348   dpe_orig = my_tri->edge_list.end();
00349 
00350   while ( dpe != dpe_orig ) {
00351     edge = *dpe;
00352     if ( edge->edge_type == BDRY_EDGE ) {
00353       dpe = my_tri->edge_list.erase(dpe);
00354       orphaned_edges.push_back(edge);
00355       dpe_orig = my_tri->edge_list.end();
00356     } else {
00357       dpe++;
00358     }
00359   }  
00360   int newv0, newv1, newv0_type, newv1_type;
00361   newv0_type = newv1_type = EDGE_0;
00362   newv0 = my_tri->v0;  
00363   dp = edge0_list.begin();
00364   while ( dp != edge0_list.end() ) {
00365     mypair = *dp;
00366     newv1 = mypair.second;
00367       //  It is possible to get the same vert more than once in the list.
00368       //  After sorting, they will be adjacent.  The following if statement
00369       //  causes duplicate verts to be used only once.
00370     if ( newv0 != newv1 ) {
00371       add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type);
00372       newv0 = newv1;
00373     }
00374     dp++;
00375   }
00376   newv1 = my_tri->v1;
00377   if ( newv0 != newv1 )
00378     add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type);
00379 
00380   newv0_type = newv1_type = EDGE_1;
00381   newv0 = my_tri->v1;  
00382   dp = edge1_list.begin();
00383   while ( dp != edge1_list.end() ) {
00384     mypair = *dp;
00385     newv1 = mypair.second;
00386     if ( newv0 != newv1 ) {
00387       add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type);
00388       newv0 = newv1;
00389     }
00390     dp++;
00391   }  
00392   newv1 = my_tri->v2;
00393   if ( newv0 != newv1 )
00394     add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type);
00395 
00396   newv0_type = newv1_type = EDGE_2;
00397   newv0 = my_tri->v2;  
00398   dp = edge2_list.begin();
00399   while ( dp != edge2_list.end() ) {
00400     mypair = *dp;
00401     newv1 = mypair.second;
00402     if ( newv0 != newv1 ) {
00403       add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type);
00404       newv0 = newv1;
00405     }
00406     dp++;
00407   }
00408   newv1 = my_tri->v0;
00409   if ( newv0 != newv1 )
00410     add_this_bdry_edge(newv0,newv1,newv0_type,newv1_type); 
00411     
00412   return status;
00413   
00414 }
00415 
00416 void FBRetriangulate::add_this_bdry_edge(int v0, int v1, int v0_type,
00417                                          int v1_type) 
00418 {
00419 //  Add an edge if it doesn't already exist.
00420   std::vector<FB_Edge*>::iterator dpe;
00421   bool ifoundit;
00422   FB_Edge *edge, *new_edge;
00423 
00424   ifoundit = false;
00425   dpe = my_tri->edge_list.begin();
00426   while ( dpe != my_tri->edge_list.end() ) {
00427     edge = *dpe;
00428     if ( ( ((int)edge->v0 == v0) && ((int)edge->v1 == v1) ) ||
00429          ( ((int)edge->v0 == v1) && ((int)edge->v1 == v0) ) ) {
00430       ifoundit = true; 
00431       break;
00432     }
00433     dpe++;
00434   }
00435   if ( ifoundit == false ) {
00436     new_edge = new FB_Edge(v0,v1,v0_type,v1_type,true);
00437     new_edge->edge_type = BDRY_EDGE;
00438     my_tri->edge_list.push_back(new_edge);  
00439   }
00440 }
00441 
00442 void FBRetriangulate::make_vert_list()
00443 {
00444   unsigned int i;
00445   int v0, v1;
00446   FB_Edge *edge;
00447   VertexStuff *vstuff;
00448   std::vector<FB_Edge*>::iterator dpe;
00449   bool foundv0, foundv1; 
00450 
00451   dpe = my_tri->edge_list.begin();
00452   while ( dpe != my_tri->edge_list.end() ) {
00453     edge = *dpe;
00454     dpe++;
00455     v0 = edge->v0;
00456     v1 = edge->v1;
00457     foundv0 = foundv1 = false;
00458     
00459     for ( i = 0; i < vertstufflist.size(); i++ ) {
00460       if ( vertstufflist[i]->v0 == v0 ) {
00461         vertstufflist[i]->edge_list.push_back(edge);
00462         foundv0 = true;
00463       }
00464       if ( vertstufflist[i]->v0 == v1 ) {
00465         vertstufflist[i]->edge_list.push_back(edge);
00466         foundv1 = true;
00467       }
00468     }
00469     if ( foundv0 == false ) {
00470       vstuff = new VertexStuff(v0, edge->v0_type,verts[v0]->coord[p_dir],
00471                                verts[v0]->coord[s_dir]);
00472       vstuff->edge_list.push_back(edge);
00473       vertstufflist.push_back(vstuff);
00474     }
00475     if ( foundv1 == false ) {
00476       vstuff = new VertexStuff(v1, edge->v1_type,verts[v1]->coord[p_dir],
00477                                verts[v1]->coord[s_dir]);
00478       vstuff->edge_list.push_back(edge);
00479       vertstufflist.push_back(vstuff);
00480     }
00481     
00482   } 
00483   
00484   std::vector<VertexStuff* >::iterator vitbegin, vitend;
00485   vitbegin = vertstufflist.begin();
00486   vitend = vertstufflist.end();
00487 
00488   std::sort(vitbegin,vitend,vertstuffcompfn_less());  
00489 
00490 }
00491 
00492 void FBRetriangulate::add_edge_up(int v0, int seq)
00493 {
00494   unsigned int i, k;
00495   int v1, v0test, v1test, v0type, v1type;
00496   bool foundit, crossed;
00497   FB_Edge *edge;
00498   std::vector<FB_Edge*>::iterator dpe;
00499 
00500   v1 = 0; //To keep compiler from warning that v1 might be used uninitialized.
00501   v1type = 0;
00502   v0type = vertstufflist[seq]->v0type;
00503   for ( k = seq+1; k < vertstufflist.size(); k++ ) {
00504     foundit = true;
00505     v1 = vertstufflist[k]->v0;
00506     if ( fabs(verts[v1]->coord[p_dir]-verts[v0]->coord[p_dir]) < EPSILON )
00507       continue;
00508     v1type = vertstufflist[k]->v0type;
00509       //  v0 to v1 is the putative new edge.  Test whether it crosses any
00510       //  existing internal edge.  Any such internal edge that it crosses
00511       //  has to have an endpoint higher than v1.  If v1 is the next-to-top
00512       //  vertex, use it.
00513     if ( k == vertstufflist.size() - 1 ) {
00514       foundit = true;
00515       add_tri_edge(v0,v1,v0type,v1type);
00516       return;      
00517     }
00518     for ( i = k+1; i < vertstufflist.size()-1; i++ ) {
00519       dpe = vertstufflist[i]->edge_list.begin();
00520 
00521       while ( dpe != vertstufflist[i]->edge_list.end() ) {
00522         edge = *dpe;
00523         dpe++;
00524         if ( (edge->v0_type != INTERIOR_VERT) &&
00525              (edge->v1_type != INTERIOR_VERT) )
00526           continue;
00527         v0test = edge->v0;
00528         v1test = edge->v1;
00529         if ( (v0test == v1) || (v1test == v1) ) continue;
00530         crossed = test_for_crossing(v0,v1,v0test,v1test);
00531         if ( crossed == true ) {
00532           foundit = false;
00533           break;
00534         }
00535       }
00536       if ( foundit == false ) break;
00537     }  //  end of "for (i ...."
00538     if ( foundit == true ) break;
00539   }
00540   add_tri_edge(v0,v1,v0type,v1type);
00541 }
00542 
00543 void FBRetriangulate::add_edge_down(int v0, int seq)
00544 {
00545   int i, k;
00546   int v1, v0test, v1test, v0type, v1type;
00547   bool foundit, crossed;
00548   FB_Edge *edge;
00549   std::vector<FB_Edge*>::iterator dpe;
00550 
00551   v1 = 0;  //To keep compiler from warning that v1 might be used uninitialized.
00552   v1type = 0;
00553   v0type = vertstufflist[seq]->v0type;
00554   for ( k = seq-1; k > -1; k-- ) {
00555     foundit = true;
00556     v1 = vertstufflist[k]->v0;
00557     if ( fabs(verts[v1]->coord[p_dir]-verts[v0]->coord[p_dir]) < EPSILON ) continue;
00558     v1type = vertstufflist[k]->v0type;
00559       //  v0 to v1 is the putative new edge.  Test whether it crosses
00560       //  any existing internal edge.  Any such internal edge that it
00561       //  crosses has to have an endpoint higher than v1.  If v1 is the
00562       //  next-to-top vertex, use it.
00563     if ( k == 0 ) {
00564       foundit = true;
00565       add_tri_edge(v0,v1,v0type,v1type);
00566       return;      
00567     }    
00568     for ( i = k-1; i > 0; i-- ) {
00569       dpe = vertstufflist[i]->edge_list.begin();
00570 
00571       while ( dpe != vertstufflist[i]->edge_list.end() ) {
00572         edge = *dpe;
00573         dpe++;
00574         if ( (edge->v0_type != INTERIOR_VERT) &&
00575              (edge->v1_type != INTERIOR_VERT) )
00576           continue;
00577         v0test = edge->v0;
00578         v1test = edge->v1;
00579         if ( (v0test == v1) || (v1test == v1) ) continue;
00580         crossed = test_for_crossing(v0,v1,v0test,v1test);
00581         if ( crossed == true ) {
00582           foundit = false;
00583           break;
00584         }
00585       }
00586       if ( foundit == false ) break;
00587     }  //  end of "for (i ...."
00588     if ( foundit == true ) break;
00589   }
00590   add_tri_edge(v0,v1,v0type,v1type);
00591 }
00592   
00593 bool FBRetriangulate::test_for_crossing(int v0, int v1, int v2, int v3)
00594 {
00595   double x0, y0, x1, y1, x2, y2, x3, y3, dxa, dya, dxb, dyb, p01x, p01y;
00596   double product, dasq, dbsq, prodsq;
00597   double s, t;
00598 
00599   x0 = verts[v0]->coord[p_dir]; y0 = verts[v0]->coord[s_dir];
00600   x1 = verts[v1]->coord[p_dir]; y1 = verts[v1]->coord[s_dir];
00601   x2 = verts[v2]->coord[p_dir]; y2 = verts[v2]->coord[s_dir];
00602   x3 = verts[v3]->coord[p_dir]; y3 = verts[v3]->coord[s_dir];
00603   dxa = x1 - x0; dya = y1 - y0;
00604   dxb = x3 - x2; dyb = y3 - y2;
00605   
00606   product = dxa*dyb - dya*dxb;
00607   dasq = dxa*dxa + dya*dya;
00608   dbsq = dxb*dxb + dyb*dyb;
00609   prodsq = product*product;
00610   
00611   if ( prodsq > EPSILON2*dasq*dbsq ) {
00612     p01x = x2 - x0;
00613     p01y = y2 - y0;
00614     s = (p01x*dyb - p01y*dxb)/product;
00615     if ( (s < 0.0) || (s > 1.0) ) return false; 
00616   
00617     t = (p01x*dya - p01y*dxa)/product;
00618     if ( (t < 0.0) || (t > 1.0) ) return false; 
00619     
00620   }
00621   
00622   return true;
00623 }
00624 
00625 void FBRetriangulate::add_tri_edge(int v0, int v1, int v0_type, int v1_type)
00626 {
00627   FB_Edge *edge;
00628   unsigned int i;
00629 
00630   edge = new FB_Edge(v0,v1,v0_type,v1_type,false);
00631   edge->edge_type = INTERIOR_EDGE;
00632   my_tri->edge_list.push_back(edge);
00633   for ( i = 0; i < vertstufflist.size(); i++ ) {
00634     if ( vertstufflist[i]->v0 == v0 ) {
00635       vertstufflist[i]->edge_list.push_back(edge);
00636     }
00637     if ( vertstufflist[i]->v0 == v1 ) {
00638       vertstufflist[i]->edge_list.push_back(edge);
00639       v1_type = vertstufflist[i]->v0type;
00640     }
00641   }
00642   
00643 }
00644 
00645 bool FBRetriangulate::get_a_chain(std::vector<int> *chainlist)
00646 {
00647   bool status;
00648   int vthis, vprev, vstart;
00649   int direction;
00650   FB_Edge *edge;
00651   std::vector<FB_Edge*>::iterator dpe;
00652   status = false;
00653   dpe = my_tri->edge_list.begin();
00654     // to keep compiler from warning that edge might be used uninitialized
00655   edge = 0;
00656   
00657   while ( dpe !=  my_tri->edge_list.end() ) {
00658     edge = *dpe;
00659     dpe++;
00660       //  Skip edges that are not interior edges.
00661 //    if ( edge->v0_type == edge->v1_type ) continue;
00662     if ( edge->edge_type == BDRY_EDGE  ) continue;
00663     if ( edge->num_times < 2 ) {
00664       if ( (edge->edge_type == BDRY_EDGE) && (edge->num_times == 1) )
00665         continue;
00666       status = true;
00667       edge->num_times++;
00668       break;
00669     }
00670   }
00671   if ( status == true ) {
00672     if ( edge->num_times == 1 ) {
00673       vthis = edge->v1;
00674       vprev = edge->v0;
00675       direction = 1;
00676     } else {
00677       vthis = edge->v0;
00678       vprev = edge->v1;
00679       direction = 1;
00680     }
00681     vstart = vprev;
00682     chainlist->push_back(vthis);
00683     
00684     while ( vthis != vstart ) {
00685       get_next_vert_in_chain(vthis,vprev,direction);
00686       chainlist->push_back(vthis);
00687     }    
00688   }
00689 
00690   if ( status == false ) { // There were no interior edges.
00691                            // So just go around the perimeter.
00692     dpe = my_tri->edge_list.begin();
00693     edge = *dpe;
00694     if ( edge->num_times != 0 ) {
00695       status = false;
00696       return status;
00697     }
00698     edge->num_times = 1;
00699     vthis = edge->v0;
00700     vprev = edge->v1;
00701     direction = 1;
00702     vstart = vprev;
00703     chainlist->push_back(vthis);
00704     
00705     while ( vthis != vstart ) {
00706       get_next_vert_in_chain(vthis,vprev,direction);
00707       chainlist->push_back(vthis);
00708     } 
00709     status = true;  
00710   }
00711   
00712   return status;
00713  
00714 }
00715 
00716 void FBRetriangulate::get_next_vert_in_chain(int& vthis,
00717                                              int& vprev, int direction)
00718 //  Right now, direction always equals 1.
00719 {
00720   unsigned int i;
00721   FB_Edge *edge;
00722   std::vector<FB_Edge*>::iterator dpe, dpebegin, dpeend;
00723 
00724   for ( i = 0; i < vertstufflist.size(); i++ ) {
00725     if ( vertstufflist[i]->v0 == (int)vthis ) break;
00726   }
00727   dpe = dpebegin = vertstufflist[i]->edge_list.begin();
00728   dpeend = vertstufflist[i]->edge_list.end();
00729   while ( dpe != dpeend ) {
00730     edge = *dpe;
00731     if ( (edge->v0 == vprev) || (edge->v1 == vprev) ) break;  
00732     dpe++;
00733   }
00734   if ( direction == 1 ) {
00735     if ( dpe == dpebegin )
00736       dpe = vertstufflist[i]->edge_list.end();
00737     dpe--;    
00738   } else {
00739     dpe++;
00740     if ( dpe == dpeend )
00741       dpe = dpebegin;
00742   }
00743   edge = *dpe;
00744   vprev = vthis;
00745   if ( vthis == edge->v1 ) {
00746       //  Swap the edge verts.  We do this so that next time we see this
00747       //  edge we will know that it was oriented to point in the direction
00748       //  of the previous loop.  Thus we will know to proceed in the
00749       //  opposite direction next time.
00750       //  (see edge->num_times test in get_a_chain().)
00751     unsigned int vtemp;
00752     vtemp = edge->v0;
00753     edge->v0 = edge->v1;
00754     edge->v1 = vtemp;
00755     vtemp = edge->v0_type;
00756     edge->v0_type = edge->v1_type;
00757     edge->v1_type = vtemp;
00758   }
00759   vthis = edge->v1;
00760   edge->num_times++;
00761   
00762 }
00763 
00764 void FBRetriangulate::sort_vertstufflist_edges()
00765 {
00766 //  For each vertex that has more than 2 edges, sort the edges in CCW order.
00767   unsigned int i;
00768   FB_Edge *edge;
00769   std::vector<FB_Edge*>::iterator dpe, dpebegin, dpeend;
00770   double x0, y0, x1, y1, slope;
00771   int v0;
00772   unsigned int quadrant;
00773 
00774   quadrant = 0;  //  Initialize so compiler won't warn.
00775   
00776   for ( i = 0; i < vertstufflist.size(); i++ ) {
00777     if ( vertstufflist[i]->edge_list.size() > 2 ) {
00778       v0 = vertstufflist[i]->v0;
00779       x0 = verts[v0]->coord[p_dir];
00780       y0 = verts[v0]->coord[s_dir];      
00781       dpe = vertstufflist[i]->edge_list.begin();
00782       while ( dpe != vertstufflist[i]->edge_list.end() ) {
00783         edge = *dpe;
00784         dpe++;
00785         if ( edge->v0 == v0 ) {
00786           x1 = verts[edge->v1]->coord[p_dir];
00787           y1 = verts[edge->v1]->coord[s_dir];
00788         } else {
00789           x1 = verts[edge->v0]->coord[p_dir];
00790           y1 = verts[edge->v0]->coord[s_dir];
00791         }
00792         if ( fabs(x1-x0) < EPSILON ) {
00793           if ( y1 > y0 ) {
00794             slope = CUBIT_DBL_MAX;
00795             quadrant = 1;
00796           } else {
00797             slope = -CUBIT_DBL_MAX;
00798             quadrant = 4;
00799           }
00800         } else {
00801           slope = (y1-y0)/(x1-x0);
00802           if ( (x1 >= x0) && (y1 >= y0) ) quadrant = 1;
00803           else if ( (x1 < x0) && (y1 > y0) ) quadrant = 2;
00804           else if ( (x1 <= x0) && (y1 <= y0) ) quadrant = 3;
00805           else if ( (x1 > x0) && (y1 < y0) ) quadrant = 4;
00806         }
00807         edge->slope = slope;
00808         edge->quadrant = quadrant;
00809       }
00810         //  Now sort the edge list by the value of quadrant or slope.
00811 //      vertstufflist[i]->edge_list.sort(edgecompfn_less());
00812 //      vertstufflist[i]->edge_list.sort(edgecf_less);
00813       dpebegin = vertstufflist[i]->edge_list.begin();
00814       dpeend = vertstufflist[i]->edge_list.end();   
00815 //      if ( winding == CW )   
00816       std::sort(dpebegin,dpeend,edgecompfn_less());
00817 //      else
00818 //        std::sort(dpebegin,dpeend,edgecompfn_more());
00819       
00820     }
00821 
00822     dpe = vertstufflist[i]->edge_list.begin();
00823     while ( dpe != vertstufflist[i]->edge_list.end() ) {
00824       edge = *dpe;
00825       dpe++;
00826     }
00827   
00828   
00829   }
00830 }
00831 
00832 void FBRetriangulate::classify_edges()
00833 {
00834 //  Flag boundary edges as such.
00835   FB_Edge *edge;
00836   std::vector<FB_Edge*>::iterator dpe;
00837   int type;
00838 
00839   dpe = my_tri->edge_list.begin();
00840   while ( dpe != my_tri->edge_list.end() ) {
00841     edge = *dpe;
00842     
00843     if ( (edge->v0_type != INTERIOR_VERT) &&
00844          (edge->v1_type != INTERIOR_VERT) ) {
00845       if ( edge->v0_type == edge->v1_type )
00846         edge->edge_type = BDRY_EDGE;
00847       else {
00848         type = UNKNOWN;
00849         switch( edge->v0_type ) {
00850           case VERTEX_0:
00851             if ( edge->v1_type != EDGE_1 ) {
00852               type = BDRY_EDGE;
00853               if ( edge->edge_type == INTERSECTION_EDGE) {
00854                 if ( edge->v1_type == VERTEX_1 ) 
00855                   my_tri->cubitedge0index = INTERSECTION_EDGE;
00856                 else if ( edge->v1_type == VERTEX_2 ) 
00857                   my_tri->cubitedge2index = INTERSECTION_EDGE;
00858               }
00859             }          
00860             break;
00861           case VERTEX_1:
00862             if ( edge->v1_type != EDGE_2 ) {
00863               type = BDRY_EDGE;
00864               if ( edge->edge_type == INTERSECTION_EDGE) {
00865                 if ( edge->v1_type == VERTEX_2 ) 
00866                   my_tri->cubitedge1index = INTERSECTION_EDGE;
00867                 else if ( edge->v1_type == VERTEX_0 ) 
00868                   my_tri->cubitedge0index = INTERSECTION_EDGE;
00869               } 
00870             }          
00871             break;
00872           case VERTEX_2:
00873             if ( edge->v1_type != EDGE_0 ) {
00874               type = BDRY_EDGE;
00875               if ( edge->edge_type == INTERSECTION_EDGE) {
00876                 if ( edge->v1_type == VERTEX_1 ) 
00877                   my_tri->cubitedge1index = INTERSECTION_EDGE;
00878                 else if ( edge->v1_type == VERTEX_0 ) 
00879                   my_tri->cubitedge2index = INTERSECTION_EDGE;
00880               }
00881             }          
00882             break;
00883           case EDGE_0:
00884             if ( (edge->v1_type == VERTEX_0) ||
00885                  (edge->v1_type == VERTEX_1) ) {
00886               type = BDRY_EDGE;
00887               if ( edge->edge_type == INTERSECTION_EDGE)
00888                 my_tri->cubitedge0index = INTERSECTION_EDGE;
00889             }
00890             break;
00891           case EDGE_1:
00892             if ( (edge->v1_type == VERTEX_1) ||
00893                  (edge->v1_type == VERTEX_2) ) {
00894               type = BDRY_EDGE;
00895               if ( edge->edge_type == INTERSECTION_EDGE)
00896                 my_tri->cubitedge1index = INTERSECTION_EDGE;
00897             }
00898             break;
00899           case EDGE_2:
00900             if ( (edge->v1_type == VERTEX_2) ||
00901                  (edge->v1_type == VERTEX_0) ) {
00902               type = BDRY_EDGE;
00903               if ( edge->edge_type == INTERSECTION_EDGE)
00904                 my_tri->cubitedge2index = INTERSECTION_EDGE;
00905             }
00906             break;
00907         }
00908         if ( type == BDRY_EDGE ) edge->edge_type = BDRY_EDGE;      
00909       }        
00910     } // else edge->edge_type = INTERIOR_EDGE;
00911     dpe++;
00912   }
00913        
00914 }
00915 
00916 void FBRetriangulate::get_edge_indices(int v0, int v1, int v2, int parent, 
00917                                        int &e0index, int &e1index,
00918                                        int &e2index)
00919 {
00920   std::vector<FB_Edge*>::iterator itt;
00921   std::vector<FB_Triangle*>::iterator itp;
00922   FB_Edge *edge;
00923   int e_v0, e_v1, e_type, e_v0type, e_v1type;
00924 
00925   itp = tris->begin();  
00926   itp += parent;
00927   itt = (*itp)->edge_list.begin();
00928   while ( itt != (*itp)->edge_list.end() ) {
00929     edge = *itt;
00930     e_v0 = edge->v0;
00931     e_v1 = edge->v1;
00932     e_v0type = edge->v0_type;
00933     e_v1type = edge->v1_type;
00934     e_type = edge->edge_type;
00935     if ( ( (e_v0 == v0) && (e_v1 == v1) ) ||
00936          ( (e_v0 == v1) && (e_v1 == v0) ) ) {
00937       if ( e_type == INTERSECTION_EDGE )
00938         e0index = INTERSECTION_EDGE; 
00939       else if ( e_type == INTERIOR_EDGE )
00940         e0index = 0;  
00941       else if ( (e_v0type == EDGE_0) || (e_v1type == EDGE_0) )
00942         e0index = (*itp)->cubitedge0index;
00943       else if ( (e_v0type == EDGE_1) || (e_v1type == EDGE_1) )
00944         e0index = (*itp)->cubitedge1index;
00945       else if ( (e_v0type == EDGE_2) || (e_v1type == EDGE_2) )
00946         e0index = (*itp)->cubitedge2index;      
00947       else if ( ( (e_v0type == VERTEX_0) || (e_v0type == VERTEX_1) ) &
00948                 ( (e_v1type == VERTEX_0) || (e_v1type == VERTEX_1) ) )
00949         e0index = (*itp)->cubitedge0index;
00950       else if ( ( (e_v0type == VERTEX_1) || (e_v0type == VERTEX_2) ) &
00951                 ( (e_v1type == VERTEX_1) || (e_v1type == VERTEX_2) ) )
00952         e1index = (*itp)->cubitedge1index;
00953       else if ( ( (e_v0type == VERTEX_2) || (e_v0type == VERTEX_0) ) &
00954                 ( (e_v1type == VERTEX_2) || (e_v1type == VERTEX_0) ) )
00955         e2index = (*itp)->cubitedge2index;                
00956     }
00957     else if ( ( (e_v0 == v1) && (e_v1 == v2) ) ||
00958               ( (e_v0 == v2) && (e_v1 == v1) ) ) {
00959       if ( e_type == INTERSECTION_EDGE )
00960         e1index = INTERSECTION_EDGE; 
00961       else if ( e_type == INTERIOR_EDGE )
00962         e1index = 0;    
00963       else if ( (e_v0type == EDGE_0) || (e_v1type == EDGE_0) )
00964         e1index = (*itp)->cubitedge0index;
00965       else if ( (e_v0type == EDGE_1) || (e_v1type == EDGE_1) )
00966         e1index = (*itp)->cubitedge1index;
00967       else if ( (e_v0type == EDGE_2) || (e_v1type == EDGE_2) )
00968         e1index = (*itp)->cubitedge2index;
00969       else if ( ( (e_v0type == VERTEX_0) || (e_v0type == VERTEX_1) ) &
00970                 ( (e_v1type == VERTEX_0) || (e_v1type == VERTEX_1) ) )
00971         e0index = (*itp)->cubitedge0index;
00972       else if ( ( (e_v0type == VERTEX_1) || (e_v0type == VERTEX_2) ) &
00973                 ( (e_v1type == VERTEX_1) || (e_v1type == VERTEX_2) ) )
00974         e1index = (*itp)->cubitedge1index;
00975       else if ( ( (e_v0type == VERTEX_2) || (e_v0type == VERTEX_0) ) &
00976                 ( (e_v1type == VERTEX_2) || (e_v1type == VERTEX_0) ) )
00977         e2index = (*itp)->cubitedge2index;
00978     }
00979     else if ( ( (e_v0 == v2) && (e_v1 == v0) ) ||
00980               ( (e_v0 == v0) && (e_v1 == v2) ) ) {
00981       if ( e_type == INTERSECTION_EDGE )
00982         e2index = INTERSECTION_EDGE; 
00983       else if ( e_type == INTERIOR_EDGE )
00984         e2index = 0;    
00985       else if ( (e_v0type == EDGE_0) || (e_v1type == EDGE_0) )
00986         e2index = (*itp)->cubitedge0index;
00987       else if ( (e_v0type == EDGE_1) || (e_v1type == EDGE_1) )
00988         e2index = (*itp)->cubitedge1index;
00989       else if ( (e_v0type == EDGE_2) || (e_v1type == EDGE_2) )
00990         e2index = (*itp)->cubitedge2index;
00991       else if ( ( (e_v0type == VERTEX_0) || (e_v0type == VERTEX_1) ) &
00992                 ( (e_v1type == VERTEX_0) || (e_v1type == VERTEX_1) ) )
00993         e0index = (*itp)->cubitedge0index;
00994       else if ( ( (e_v0type == VERTEX_1) || (e_v0type == VERTEX_2) ) &
00995                 ( (e_v1type == VERTEX_1) || (e_v1type == VERTEX_2) ) )
00996         e1index = (*itp)->cubitedge1index;
00997       else if ( ( (e_v0type == VERTEX_2) || (e_v0type == VERTEX_0) ) &
00998                 ( (e_v1type == VERTEX_2) || (e_v1type == VERTEX_0) ) )
00999         e2index = (*itp)->cubitedge2index;
01000     }
01001     
01002     itt++;
01003   }
01004 
01005 }
01006 
01007                                        
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines