cgma
FBTiler.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 <vector>
00030 #include "CubitMessage.hpp"
00031 #include "FBTiler.hpp"
00032 
00033 FBTiler::FBTiler(std::vector<FB_Coord *>& my_verts, int pd, int sd, int sequence,
00034                  double a, double b, double c,
00035                  std::vector<int> *tri_list)
00036 {
00037 
00038   verts = my_verts;
00039   p_dir = pd;
00040   s_dir = sd;
00041   xnorm = a;
00042   ynorm = b;
00043   znorm = c;
00044   parent = sequence;
00045   my_tri_list = tri_list;
00046 }
00047 
00048 FBTiler::~FBTiler()
00049 {
00050 
00051 }
00052 
00053 CubitStatus FBTiler::Tile_Poly(std::vector<int> *coordlist)
00054 {
00055 int v0, v1, v2, v3;
00056 std::vector<int>::iterator it, itmin, itmax;
00057 double min_p, min_s, max_p, max_s, test_p, test_s;
00058 CubitStatus status;
00059 
00060   status = CUBIT_SUCCESS;
00061   
00062   it = coordlist->begin();
00063   
00064   if ( coordlist->size() == 3 ) {
00065     v1 = *it++;  v2 = *it++; v3 = *it++;  
00066     add_triangle(v1, v2, v3);
00067     return status;
00068   }
00069 
00070   //  Get location in list of min and max coord in p_dir.
00071   min_p = min_s = CUBIT_DBL_MAX;
00072   max_p = max_s = -CUBIT_DBL_MAX;
00073   while ( it != coordlist->end() ) {
00074     v0 = *it;
00075     test_p = verts[v0]->coord[p_dir];
00076     test_s = verts[v0]->coord[s_dir];
00077 //    if ( test_p < min_p ) {
00078     if ( (test_p - min_p) <= -EPSILON ) {
00079       itmin = it;
00080       min_p = test_p;
00081       if ( test_s < min_s ) 
00082         min_s = test_s;
00083     } 
00084 //    if ( test_p > max_p ) {
00085     if ( (test_p - max_p) >= EPSILON ) {
00086       itmax = it;
00087       max_p = test_p;
00088       if ( test_s > max_s )
00089         max_s = test_s;      
00090     }
00091     if ( fabs(test_p - min_p) < EPSILON ) {
00092       if ( test_s < min_s ) {
00093         itmin = it;
00094         min_s = test_s;
00095       }
00096     } 
00097     if ( fabs(test_p - max_p) < EPSILON ) {
00098       if ( test_s > max_s ) {
00099         itmax = it;
00100         max_s = test_s;
00101       }    
00102     }
00103 
00104     it++;
00105   }  
00106 
00107   //  Now we've got to start at the min coord and make a single sorted list  
00108   //  of the coords, from min to max, noting whether the coord is on the
00109   //  LEFT or RIGHT chain, or -- for the top and bottom coords -- on BOTH chains.
00110 
00111   sortedchainlist.clear();  
00112   
00113   std::vector<int>::iterator itleft, itright, itlistbegin, itlistend;
00114   it = itleft = itright = itmin;  
00115   itlistbegin = coordlist->begin();
00116   itlistend = coordlist->end();
00117   
00118   decrement_list_ptr(itleft,itlistbegin,itlistend);
00119   increment_list_ptr(itright,itlistbegin,itlistend); 
00120   
00121   v0 = *itmin;
00122   int v0left, v0right;
00123   double leftpdircoord, rightpdircoord;
00124   FBTilerChainVert *cvert;
00125   cvert = new FBTilerChainVert(v0,BOTH);
00126   sortedchainlist.push_back(cvert);
00127     v0left = *itleft; v0right = *itright;  
00128     leftpdircoord = verts[v0left]->coord[p_dir];
00129     rightpdircoord = verts[v0right]->coord[p_dir];
00130 
00131   while(1) {  
00132     while ( (leftpdircoord - rightpdircoord) > -EPSILON2 ) {
00133       if ( itright == itmax ) break;
00134       v0 = v0right;
00135       cvert = new FBTilerChainVert(v0,LEFT);
00136       sortedchainlist.push_back(cvert);
00137       increment_list_ptr(itright,itlistbegin,itlistend); 
00138       v0right = *itright;
00139       rightpdircoord = verts[v0right]->coord[p_dir];
00140 
00141     }
00142     while ( (rightpdircoord - leftpdircoord) > -EPSILON2 ) {
00143       if ( itleft == itmax ) break;
00144       v0 = v0left;
00145       cvert = new FBTilerChainVert(v0,RIGHT);
00146       sortedchainlist.push_back(cvert);
00147       decrement_list_ptr(itleft,itlistbegin,itlistend);
00148       v0left = *itleft;
00149       leftpdircoord = verts[v0left]->coord[p_dir];
00150 
00151     }    
00152     if ( (itright == itmax) || (itleft == itmax) ) break;
00153   }
00154 
00155 //  In case we have some horizontal edges (wrt p_dir) left on the list,
00156 //  the next two while statements will get them.
00157   while ( itright != itmax ) {
00158       v0 = v0right;
00159       cvert = new FBTilerChainVert(v0,LEFT);
00160       sortedchainlist.push_back(cvert);
00161       increment_list_ptr(itright,itlistbegin,itlistend); 
00162       v0right = *itright;
00163   } 
00164   while ( itleft != itmax ) {
00165        v0 = v0left;
00166       cvert = new FBTilerChainVert(v0,RIGHT);
00167       sortedchainlist.push_back(cvert);
00168       decrement_list_ptr(itleft,itlistbegin,itlistend);
00169       v0left = *itleft;
00170   } 
00171 
00172 //  Put the top coord on the list.
00173 
00174   v0 = *itmax;
00175   cvert = new FBTilerChainVert(v0,BOTH);
00176   sortedchainlist.push_back(cvert);
00177   status = retriangulate(coordlist);
00178 std::vector<FBTilerChainVert*>::iterator itv, itvlast;  
00179 
00180   itvlast = sortedchainlist.end();
00181   itv = sortedchainlist.begin();
00182   while ( itv != itvlast ) {
00183     delete *itv;   
00184     itv++;
00185   }
00186   sortedchainlist.clear();
00187   return status;
00188 }
00189 
00190 int FBTiler::add_triangle(int v1, int v2, int v3)
00191 {
00192 int status;
00193 
00194   status = 1;
00195 
00196   //  See if the winding is OK by comparing the normal to the parent's normal.
00197   double xn, yn, zn;
00198   double x1, y1, z1, x2, y2, z2, x3, y3, z3;
00199   double ux, uy, uz, vx, vy, vz, product;
00200 
00201   x1 = verts[v1]->coord[0]; y1 = verts[v1]->coord[1]; z1 = verts[v1]->coord[2];
00202   x2 = verts[v2]->coord[0]; y2 = verts[v2]->coord[1]; z2 = verts[v2]->coord[2];
00203   x3 = verts[v3]->coord[0]; y3 = verts[v3]->coord[1]; z3 = verts[v3]->coord[2];
00204   ux = x2 - x1; uy = y2 - y1; uz = z2 - z1;
00205   vx = x3 - x1; vy = y3 - y1; vz = z3 - z1;
00206   xn = uy*vz - uz*vy; yn = uz*vx - ux*vz; zn = ux*vy - uy*vx;
00207   product = xnorm*xn + ynorm*yn + znorm*zn;
00208   if ( product < 0.0 ) {  //  Must reverse the connections.
00209     int itemp;
00210     itemp = v1; v1 = v2; v2 = itemp;
00211   }
00212   my_tri_list->push_back(v1);
00213   my_tri_list->push_back(v2);
00214   my_tri_list->push_back(v3);
00215  
00216   return status;
00217 
00218 }
00219 
00220 void FBTiler::dud_from_coord_list(int val, std::vector<int> *ivec)
00221 {
00222 std::vector<int>::iterator itv;
00223 bool ifoundit;
00224 //  Also decrements ivec->size() if val was found.
00225 
00226   ifoundit = false;
00227   itv = ivec->begin();
00228   do {
00229     if ( ifoundit == false ) {
00230       if ( *itv == val ) {
00231         ifoundit = true;
00232         break;  
00233       }
00234     }
00235     itv++;
00236   } while ( itv != ivec->end() );
00237   if ( ifoundit == true ) {
00238     ivec->erase(itv);
00239   } 
00240 
00241 }
00242 
00243 bool FBTiler::reflex_angle(int v0, int v1, int v2, std::vector<int> *coordlist)
00244 {
00245 double v0x, v0y, v1x, v1y, v2x, v2y, xbary, ybary;
00246 
00247   v0x = verts[v0]->coord[s_dir]; 
00248   v0y = verts[v0]->coord[p_dir]; 
00249   v1x = verts[v1]->coord[s_dir]; 
00250   v1y = verts[v1]->coord[p_dir]; 
00251   v2x = verts[v2]->coord[s_dir]; 
00252   v2y = verts[v2]->coord[p_dir]; 
00253 //  Are the three points colinear?
00254 
00255   if ( fabs((v1x-v0x)*(v2y-v0y) - (v2x-v0x)*(v1y-v0y)) < EPSILON2 ) 
00256     return true;
00257         
00258 //  Get coordinates of barycenter of the triangle formed by 
00259 //  v0, v1, and v2
00260 
00261   xbary = (v0x + v1x + v2x)/3.;
00262   ybary = (v0y + v1y + v2y)/3.;
00263   
00264 //  Do a point-in-polygon test on this point and the polygon formed
00265 //  by the chain.  If point is in polygon, angle is not reflex; 
00266 //  otherwise, angle is reflex.  Point in polygon algorithm from
00267 //  Schneider and Eberly, Geometric Tools for COmputer Graphics, 
00268 //  Sec. 13.3.
00269 
00270 bool inside = false;
00271 std::vector<int>::iterator itv;
00272 int u0, u1;
00273 
00274   itv = coordlist->begin();
00275   while ( itv != coordlist->end() ) {
00276     u0 = *itv++;
00277     if ( itv == coordlist->end() )
00278       u1 = *(coordlist->begin());
00279     else 
00280       u1 = *itv;
00281     if ( ybary < verts[u1]->coord[p_dir] ) {
00282     //  u1 is above the hroizontal ray in s_dir from the barycenter.
00283       if ( verts[u0]->coord[p_dir] <= ybary ) {
00284       //  u0 is on or below the ray.
00285         if ( ( (ybary - verts[u0]->coord[p_dir])*
00286            (verts[u1]->coord[s_dir] -verts[ u0]->coord[s_dir]) ) >
00287          ( (xbary - verts[u0]->coord[s_dir])*
00288            (verts[u1]->coord[p_dir] - verts[u0]->coord[p_dir]) ) )
00289       inside = !inside; 
00290     }    
00291     } else if ( ybary < verts[u0]->coord[p_dir] ) {
00292     //  U1 is on or below the ray; u0 is above the ray.
00293         if ( ( (ybary - verts[u0]->coord[p_dir])*
00294            (verts[u1]->coord[s_dir] - verts[u0]->coord[s_dir]) ) <
00295          ( (xbary - verts[u0]->coord[s_dir])*
00296            (verts[u1]->coord[p_dir] - verts[u0]->coord[p_dir]) ) )
00297       inside = !inside; 
00298     }   
00299   }
00300 
00301   return !inside;  
00302   
00303 }
00304 
00305 CubitStatus FBTiler::retriangulate(std::vector<int> *coordlist)
00306 {
00307 int stacktopadjacency,  this_adjacency, adjacency_case;
00308 std::vector<FBTilerChainVert*> vstack;
00309 std::vector<FBTilerChainVert*>::iterator itv, itvlast;
00310 int v0, v1, v2;
00311 unsigned int i, numtris, totaltris;
00312 
00313 //  Put the first two verts on vstack.
00314   itv = sortedchainlist.begin();
00315   vstack.push_back(*itv++);  
00316   vstack.push_back(*itv);
00317   stacktopadjacency = (*itv)->whichchain;
00318   itvlast = itv;
00319   numtris = 0;
00320   totaltris = sortedchainlist.size() - 2;
00321   
00322   while(1) {
00323     
00324     if ( numtris == totaltris ) break;
00325     *itv++;
00326     if ( itv == sortedchainlist.end() ) {
00327       PRINT_ERROR("Tiler error in FacetBoolean:  premature end of list.\n");
00328       return CUBIT_FAILURE;
00329     }
00330  
00331     this_adjacency = (*itv)->whichchain;
00332     v0 = (*itv)->v0;
00333 
00334     adjacency_case = get_adjacency(stacktopadjacency,this_adjacency);
00335 
00336     switch ( adjacency_case ) {
00337   
00338       case 1:
00339       {
00340         v1 = vstack[0]->v0;
00341         for ( i = 1; i < vstack.size(); i++ ) {
00342           v2 = vstack[i]->v0;
00343           add_triangle(v0,v1,v2);
00344           dud_from_coord_list(v1,coordlist);
00345           numtris += 1;
00346           v1 = v2;
00347         }   
00348         vstack.clear();
00349         vstack.push_back(*itvlast);
00350         vstack.push_back(*itv);
00351         stacktopadjacency = (*itv)->whichchain;
00352         itvlast = itv;
00353 
00354         break;
00355       }
00356 
00357       case 2:
00358       {
00359         while ( vstack.size() > 1 ) {
00360       int size = vstack.size() - 1;
00361       v1 = vstack[size]->v0;
00362 //      int v1chain = vstack[size]->whichchain;
00363       v2 = vstack[size-1]->v0;
00364 //      bool isreflex = reflex_angle(v0,v1,v2,v1chain);
00365       bool isreflex = reflex_angle(v0,v1,v2,coordlist);
00366       if ( isreflex == true ) break;
00367       add_triangle(v0,v1,v2);
00368       dud_from_coord_list(v1,coordlist);
00369       numtris += 1;
00370       vstack.pop_back();
00371 
00372    }
00373     
00374    vstack.push_back(*itv);
00375    stacktopadjacency = (*itv)->whichchain;
00376         itvlast = itv;
00377     
00378     
00379         break;
00380       }  
00381   
00382     }
00383  
00384   } // endwhile
00385 
00386   return CUBIT_SUCCESS;
00387   
00388 }
00389 
00390 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines