cgma
FBIntersect.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 <vector>
00028 #include "FBIntersect.hpp"
00029 #include "FBPolyhedron.hpp"
00030 #include "FBRetriangulate.hpp"
00031 #include "CubitMessage.hpp"
00032 #include "GfxDebug.hpp"
00033 #include "GeometryDefines.h"
00034 
00035 FBIntersect::FBIntersect()
00036 {
00037   poly1 = poly2 = 0;
00038   do_classify = false;
00039   do_edges_only = false;
00040   do_imprint = false;
00041   classify1 = classify2 = 0;
00042   body1_is_plane = body2_is_plane = false;
00043   f_c_indices1 = f_c_indices2 = 0;
00044   nothing_intersected = false;
00045 }
00046 
00047 FBIntersect::~FBIntersect()
00048 {
00049   if ( poly1 ) delete poly1;
00050   if ( poly2 ) delete poly2;
00051   if ( classify1 ) delete classify1;
00052   if ( classify2 ) delete classify2;
00053 }
00054 
00055 CubitStatus FBIntersect::intersect(const std::vector<double>& Ticoords,
00056                      const std::vector<int>& Ticonnections,
00057                      const std::vector<double>& Tjcoords,
00058                      const std::vector<int>& Tjconnections,
00059                      std::vector<int>& duddedTiFacets, 
00060                      std::vector<int>& duddedTjFacets,
00061                      std::vector<int>& newTiFacets, 
00062                      std::vector<int>& newTjFacets,
00063                      std::vector<int>& newTiFacetsIndex,
00064                      std::vector<int>& newTjFacetsIndex,
00065                      std::vector<double>& newTiPoints, 
00066                      std::vector<double>& newTjPoints,
00067                      std::vector<int>& edgesTi, 
00068                      std::vector<int>& edgesTj
00069                      )
00070 {
00071 CubitStatus status;
00072 bool boxes_intersect;
00073 
00074   status = CUBIT_SUCCESS;   
00075 
00076   if ( Ticoords.size()%3 != 0 ) {
00077     PRINT_ERROR("Bad coordinates for first part fed to FBIntersect.\n");
00078     return CUBIT_FAILURE;
00079   }
00080   if ( Tjcoords.size()%3 != 0 ) {
00081     PRINT_ERROR("Bad coordinates for second part fed to FBIntersect.\n");
00082     return CUBIT_FAILURE;
00083   }
00084   if ( Ticonnections.size()%3 != 0 ) {
00085     PRINT_ERROR("Bad connection list for first part fed to FBIntersect.\n");
00086     return CUBIT_FAILURE;
00087   }
00088   if ( Tjconnections.size()%3 != 0 ) {
00089     PRINT_ERROR("Bad connection list for second part fed to FBIntersect.\n");
00090     return CUBIT_FAILURE;
00091   }    
00092      
00093   poly1 = new FBPolyhedron;
00094   status = poly1->makepoly(Ticoords,Ticonnections,f_c_indices1);
00095   if ( status != CUBIT_SUCCESS )
00096   {
00097       return status;
00098   }
00099 
00100   poly2 = new FBPolyhedron;
00101   status = poly2->makepoly(Tjcoords,Tjconnections,f_c_indices2);
00102   if ( status != CUBIT_SUCCESS )
00103   {
00104       return status;
00105   }
00106   
00107   boxes_intersect = poly1->boxintersection(poly2);
00108   if ( boxes_intersect == false ) {
00109     nothing_intersected = true;
00110     return status;
00111   }
00112   
00113   status = pair_intersect();  
00114   if ( status != CUBIT_SUCCESS )
00115   {
00116       return status;
00117   }
00118 
00119   poly1->putnewpoints(newTiPoints);
00120   poly2->putnewpoints(newTjPoints);
00121 
00122   poly1->putedges(edgesTi);
00123   poly2->putedges(edgesTj);
00124   
00125   status = poly1->retriangulate(newTiFacets, newTiFacetsIndex);
00126   if ( status != CUBIT_SUCCESS )
00127   {
00128       return status;
00129   }
00130 
00131   status = poly2->retriangulate(newTjFacets, newTjFacetsIndex);
00132   if ( status != CUBIT_SUCCESS )
00133   {
00134       return status;
00135   }
00136 
00137   unsigned int i;
00138   for ( i = 0; i < poly1->tris.size(); i++ ) {
00139     if ( poly1->tris[i]->dudded == true ) duddedTiFacets.push_back(i);
00140   }
00141 
00142   for ( i = 0; i < poly2->tris.size(); i++ ) {
00143     if ( poly2->tris[i]->dudded == true ) duddedTjFacets.push_back(i);
00144   }
00145   
00146 //  newxxFacets is updated in FBRetriangulate::make_this_tri() after
00147 //  a new triangle has been made.  Also newxxFacetsIndex is updated if
00148 //  any new facets were made.  But this update points to the start of
00149 //  newxxFacets.  Therefore we have to add on the last member to 
00150 //  newxxFacetsIndex here. 
00151   newTiFacetsIndex.push_back(newTiFacets.size());
00152   newTjFacetsIndex.push_back(newTjFacets.size());    
00153      
00154 //  If classification of the intersected surface parts is needed, we need to
00155 //  do some things.
00156   bool other_is_planar;
00157   if ( do_classify == true ) {
00158     poly1->add_new_triangle_data();
00159     poly2->add_new_triangle_data();
00160     poly1->removeduddedtriangles();
00161     poly2->removeduddedtriangles();    
00162     classify1 = new FBClassify;
00163     classify1->SetPoly(poly1, poly2);    
00164     status = classify1->Group(1);
00165     if ( status != CUBIT_SUCCESS )
00166     {
00167       return status;
00168     }
00169     other_is_planar = body2_is_plane;
00170     status = classify1->CharacterizeGroups(1,other_is_planar);
00171     if ( status != CUBIT_SUCCESS )
00172     {
00173       return status;
00174     }
00175     classify2 = new FBClassify;
00176     classify2->SetPoly(poly1, poly2);
00177     status = classify2->Group(2);
00178     if ( status != CUBIT_SUCCESS )
00179     {
00180       return status;
00181     }
00182     other_is_planar = body1_is_plane;
00183     status = classify2->CharacterizeGroups(2,other_is_planar);
00184     if ( status != CUBIT_SUCCESS )
00185     {
00186       return status;
00187     }
00188   }  
00189   
00190   return status;
00191 
00192 }
00193 
00194 CubitStatus FBIntersect::intersect(const std::vector<double>& Ticoords,
00195                      const std::vector<int>& Ticonnections,
00196                      const std::vector<double>& Tjcoords,
00197                      const std::vector<int>& Tjconnections,
00198                      std::vector<int>& newTiFacets, 
00199                      std::vector<int>& newTjFacets, 
00200                      std::vector<int> *indices1,
00201                      std::vector<int> *indices2
00202                      )
00203 {
00204 CubitStatus status;
00205 bool boxes_intersect;
00206 
00207   status = CUBIT_SUCCESS;   
00208 
00209   if ( Ticoords.size()%3 != 0 ) {
00210     PRINT_ERROR("Bad coordinates for first part fed to FBIntersect.\n");
00211     return CUBIT_FAILURE;
00212   }
00213   if ( Tjcoords.size()%3 != 0 ) {
00214     PRINT_ERROR("Bad coordinates for second part fed to FBIntersect.\n");
00215     return CUBIT_FAILURE;
00216   }
00217   if ( Ticonnections.size()%3 != 0 ) {
00218     PRINT_ERROR("Bad connection list for first part fed to FBIntersect.\n");
00219     return CUBIT_FAILURE;
00220   }
00221   if ( Tjconnections.size()%3 != 0 ) {
00222     PRINT_ERROR("Bad connection list for second part fed to FBIntersect.\n");
00223     return CUBIT_FAILURE;
00224   }    
00225   
00226   f_c_indices1 = indices1;
00227   f_c_indices2 = indices2;
00228    
00229   poly1 = new FBPolyhedron;
00230   status = poly1->makepoly(Ticoords,Ticonnections,f_c_indices1);
00231   if ( status != CUBIT_SUCCESS )
00232   {
00233       return status;
00234   }
00235 
00236   poly2 = new FBPolyhedron;
00237   status = poly2->makepoly(Tjcoords,Tjconnections,f_c_indices2);
00238   if ( status != CUBIT_SUCCESS )
00239   {
00240       return status;
00241   }
00242 
00243   double min_angle=0.0, max_angle=0.0;
00244   int mydebug = 0;
00245   if(mydebug){
00246     poly1->min_max_angles_in_polyhedron(min_angle, max_angle);
00247     
00248     PRINT_INFO("(0) Min angle in poly1 %f, max angle  %f\n", min_angle,
00249                max_angle);
00250     poly2->min_max_angles_in_polyhedron(min_angle, max_angle);
00251     PRINT_INFO("(0) Min angle in poly2 %f, max angle  %f\n", min_angle,
00252                max_angle);
00253   }
00254   boxes_intersect = poly1->boxintersection(poly2);
00255   if ( boxes_intersect == false ) {
00256     nothing_intersected = true;
00257       //don't return because in the case of union we need to unite the
00258       //bodies even if they do not intersect
00259       //return status;
00260   }
00261   else{
00262     status = pair_intersect();
00263     if ( status != CUBIT_SUCCESS )
00264     {
00265         return status;
00266     }
00267   }
00268   
00269   if(mydebug){
00270     poly1->min_max_angles_in_polyhedron(min_angle, max_angle);
00271     PRINT_INFO("(1) Min angle in poly1 %f, max angle  %f\n", min_angle,
00272                max_angle);
00273     poly2->min_max_angles_in_polyhedron(min_angle, max_angle);
00274     PRINT_INFO("(1) Min angle in poly2 %f, max angle  %f\n", min_angle,
00275                max_angle);
00276   }
00277   
00278   status = poly1->retriangulate(newTiFacets);
00279   if ( status != CUBIT_SUCCESS )
00280   {
00281       return status;
00282   }
00283   status = poly2->retriangulate(newTjFacets);
00284   if ( status != CUBIT_SUCCESS )
00285   {
00286       return status;
00287   }
00288 
00289   if(mydebug){
00290     GfxDebug::clear();
00291     poly1->debug_draw_boundary_edges(CUBIT_BLUE_INDEX);
00292     GfxDebug::mouse_xforms();
00293     GfxDebug::clear();
00294     poly2->debug_draw_boundary_edges(CUBIT_RED_INDEX);
00295     GfxDebug::mouse_xforms();
00296     poly1->min_max_angles_in_polyhedron(min_angle, max_angle);
00297     
00298     PRINT_INFO("(2) Min angle in poly1 %f, max angle  %f\n", min_angle,
00299                max_angle);
00300     poly2->min_max_angles_in_polyhedron(min_angle, max_angle);
00301     PRINT_INFO("(2) Min angle in poly2 %f, max angle  %f\n", min_angle,
00302                max_angle);
00303   }
00304 //  If classification of the intersected surface parts is needed, we need to
00305 //  do some things.
00306   bool other_is_planar;
00307   if ( do_classify == true ) {
00308     poly1->add_new_triangle_data();
00309     poly2->add_new_triangle_data();
00310     poly1->removeduddedtriangles();
00311     poly2->removeduddedtriangles();    
00312     classify1 = new FBClassify;
00313     classify1->SetPoly(poly1, poly2);    
00314     status = classify1->Group(1);
00315     if ( status != CUBIT_SUCCESS )
00316     {
00317       return status;
00318     }
00319     other_is_planar = body2_is_plane;
00320     status = classify1->CharacterizeGroups(1,other_is_planar);
00321     if ( status != CUBIT_SUCCESS )
00322     {
00323       return status;
00324     }
00325     classify2 = new FBClassify;
00326     classify2->SetPoly(poly1, poly2);
00327     status = classify2->Group(2);
00328     if ( status != CUBIT_SUCCESS )
00329     {
00330       return status;
00331     }
00332     other_is_planar = body1_is_plane;
00333     status = classify2->CharacterizeGroups(2,other_is_planar);
00334     if ( status != CUBIT_SUCCESS )
00335     {
00336       return status;
00337     }
00338   }  
00339   
00340   return status;
00341 
00342 }
00343 
00344 CubitStatus FBIntersect::pair_intersect()
00345 {
00346 CubitStatus status;
00347 unsigned int i, j, k;
00348 int numboxesfound, *boxlist;
00349 double dtemp;
00350 
00351   boxlist = new int[poly2->tris.size()];
00352   
00353   status = CUBIT_SUCCESS;
00354   for ( i = 0; i < poly1->tris.size(); i++ ) {
00355      FB_Triangle *tri1 = poly1->tris[i];
00356      if ( (tri1->boundingbox.xmax < poly2->polyxmin) || 
00357           (tri1->boundingbox.xmin > poly2->polyxmax) ||
00358           (tri1->boundingbox.ymax < poly2->polyymin) || 
00359           (tri1->boundingbox.ymin > poly2->polyymax) ||
00360           (tri1->boundingbox.zmax < poly2->polyzmin) || 
00361           (tri1->boundingbox.zmin > poly2->polyzmax) ) continue;
00362      poly2->kdtree->box_kdtree_intersect(tri1->boundingbox,&numboxesfound,boxlist);
00363 
00364      for ( j = 0; j < (unsigned int)numboxesfound; j++ ) { 
00365         FB_Triangle *tri2 = poly2->tris[boxlist[j]];     
00366         if ( (tri1->boundingbox.xmax < tri2->boundingbox.xmin) || 
00367              (tri1->boundingbox.xmin > tri2->boundingbox.xmax) ||
00368              (tri1->boundingbox.ymax < tri2->boundingbox.ymin) || 
00369              (tri1->boundingbox.ymin > tri2->boundingbox.ymax) ||
00370              (tri1->boundingbox.zmax < tri2->boundingbox.zmin ) || 
00371              (tri1->boundingbox.zmin > tri2->boundingbox.zmax) ) continue;
00372 
00373      int mydebug = 0;
00374      if ( mydebug )
00375      {
00376          // Draw the 2 facets which are intersecting.
00377          GfxDebug::clear();
00378          poly1->debug_draw_boundary_edges(CUBIT_YELLOW_INDEX);
00379          poly2->debug_draw_boundary_edges(CUBIT_PINK_INDEX);
00380          poly1->debug_draw_fb_triangle(tri1);
00381          poly2->debug_draw_fb_triangle(tri2);
00382          GfxDebug::mouse_xforms();
00383      }
00384 
00385 //  Find the line of intersection of the two triangle planes. 
00386         linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00387         linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00388         linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00389 
00390         if ( (fabs(linecoeff[0]) < EPSILON) && 
00391              (fabs(linecoeff[1]) < EPSILON) &&  
00392              (fabs(linecoeff[2]) < EPSILON) ) {
00393             if ( do_imprint == true ) continue;
00394             //  Don't intersect nearly-coplanar triangles.
00395 
00396             // calculate the distance between the triangles.  Just because they are coplanar
00397             // does not mean we have to intersect.  If the distance between them is larger than
00398             // GEOMETRY_RESABS, then just continue on.
00399             double dist = poly1->verts[tri1->v2]->coord[0]*tri2->a +
00400                           poly1->verts[tri1->v2]->coord[1]*tri2->b +
00401                           poly1->verts[tri1->v2]->coord[2]*tri2->c + tri2->d;
00402             if ( dist > GEOMETRY_RESABS )
00403             {
00404                 continue;
00405             }
00406 
00407             //  coplanar triangles; tilt each vertex of the triangle up successively
00408             //  and then do the usual tri-tri intersection.
00409             double ta, tb, tc, td;
00410             double tx1[3],ty1[3],tz1[3];
00411 
00412             //  triangle 1
00413             ta = tri1->a; tb = tri1->b; tc = tri1->c; td = tri1->d;
00414             for ( k = 0; k < 3; k++ ) {
00415               tx1[k] = poly1->verts[tri1->v0]->coord[k];
00416               ty1[k] = poly1->verts[tri1->v1]->coord[k];
00417               tz1[k] = poly1->verts[tri1->v2]->coord[k];
00418             } 
00419 
00420             poly1->verts[tri1->v2]->coord[0] += ta;
00421             poly1->verts[tri1->v2]->coord[1] += tb;
00422             poly1->verts[tri1->v2]->coord[2] += tc;
00423             //  Compute new normal and linecoeff
00424             newplanecoefficients(poly1, tri1);
00425             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00426             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00427             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00428             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00429                          linecoeff[1]*linecoeff[1] +
00430                          linecoeff[2]*linecoeff[2]);
00431             linecoeff[0] /= dtemp;
00432             linecoeff[1] /= dtemp;
00433             linecoeff[2] /= dtemp;
00434 
00435             status = tri_tri_intersect(tri1,tri2);
00436             if ( status != CUBIT_SUCCESS )
00437             {
00438                 return status;
00439             }
00440 
00441             //  Restore values.
00442             poly1->verts[tri1->v2]->coord[0] = tz1[0];
00443             poly1->verts[tri1->v2]->coord[1] = tz1[1];
00444             poly1->verts[tri1->v2]->coord[2] = tz1[2];
00445             tri1->a = ta; tri1->b = tb; tri1->c = tc; tri1->d = td;
00446 
00447             poly1->verts[tri1->v1]->coord[0] += ta;
00448             poly1->verts[tri1->v1]->coord[1] += tb;
00449             poly1->verts[tri1->v1]->coord[2] += tc;
00450             //  Compute new normal and linecoeff
00451             newplanecoefficients(poly1, tri1);
00452             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00453             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00454             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00455             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00456                          linecoeff[1]*linecoeff[1] +
00457                          linecoeff[2]*linecoeff[2]);
00458             linecoeff[0] /= dtemp;
00459             linecoeff[1] /= dtemp;
00460             linecoeff[2] /= dtemp;
00461 
00462             status = tri_tri_intersect(tri1,tri2);
00463             if ( status != CUBIT_SUCCESS )
00464             {
00465                 return status;
00466             }
00467 
00468             //  Restore values.
00469             poly1->verts[tri1->v1]->coord[0] = ty1[0];
00470             poly1->verts[tri1->v1]->coord[1] = ty1[1];
00471             poly1->verts[tri1->v1]->coord[2] = ty1[2];
00472             tri1->a = ta; tri1->b = tb; tri1->c = tc; tri1->d = td;
00473 
00474             poly1->verts[tri1->v0]->coord[0] += ta;
00475             poly1->verts[tri1->v0]->coord[1] += tb;
00476             poly1->verts[tri1->v0]->coord[2] += tc;
00477             //  Compute new normal and linecoeff
00478             newplanecoefficients(poly1, tri1);
00479             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00480             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00481             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00482             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00483                          linecoeff[1]*linecoeff[1] +
00484                          linecoeff[2]*linecoeff[2]);
00485             linecoeff[0] /= dtemp;
00486             linecoeff[1] /= dtemp;
00487             linecoeff[2] /= dtemp;
00488 
00489             status = tri_tri_intersect(tri1,tri2);
00490             if ( status != CUBIT_SUCCESS )
00491             {
00492                 return status;
00493             }
00494 
00495             //  Restore values.
00496             poly1->verts[tri1->v0]->coord[0] = tx1[0];
00497             poly1->verts[tri1->v0]->coord[1] = tx1[1];
00498             poly1->verts[tri1->v0]->coord[2] = tx1[2];
00499             tri1->a = ta; tri1->b = tb; tri1->c = tc; tri1->d = td;
00500 
00501             //  triangle 2
00502             ta = tri2->a; tb = tri2->b; tc = tri2->c; td = tri2->d;
00503             for ( k = 0; k < 3; k++ ) {
00504               tx1[k] = poly2->verts[tri2->v0]->coord[k];
00505               ty1[k] = poly2->verts[tri2->v1]->coord[k];
00506               tz1[k] = poly2->verts[tri2->v2]->coord[k];
00507             } 
00508 
00509             poly2->verts[tri2->v2]->coord[0] += ta;
00510             poly2->verts[tri2->v2]->coord[1] += tb;
00511             poly2->verts[tri2->v2]->coord[2] += tc;
00512             //  Compute new normal and linecoeff
00513             newplanecoefficients(poly2, tri2);
00514             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00515             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00516             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00517             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00518                          linecoeff[1]*linecoeff[1] +
00519                          linecoeff[2]*linecoeff[2]);
00520             linecoeff[0] /= dtemp;
00521             linecoeff[1] /= dtemp;
00522             linecoeff[2] /= dtemp;
00523 
00524             status = tri_tri_intersect(tri1,tri2);
00525             if ( status != CUBIT_SUCCESS )
00526             {
00527                 return status;
00528             }
00529 
00530             //  Restore values.
00531             poly2->verts[tri2->v2]->coord[0] = tz1[0];
00532             poly2->verts[tri2->v2]->coord[1] = tz1[1];
00533             poly2->verts[tri2->v2]->coord[2] = tz1[2];
00534             tri2->a = ta; tri2->b = tb; tri2->c = tc; tri2->d = td;
00535 
00536             poly2->verts[tri2->v1]->coord[0] += ta;
00537             poly2->verts[tri2->v1]->coord[1] += tb;
00538             poly2->verts[tri2->v1]->coord[2] += tc;
00539             //  Compute new normal and linecoeff
00540             newplanecoefficients(poly2, tri2);
00541             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00542             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00543             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00544             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00545                          linecoeff[1]*linecoeff[1] +
00546                          linecoeff[2]*linecoeff[2]);
00547             linecoeff[0] /= dtemp;
00548             linecoeff[1] /= dtemp;
00549             linecoeff[2] /= dtemp;
00550 
00551             status = tri_tri_intersect(tri1,tri2);
00552             if ( status != CUBIT_SUCCESS )
00553             {
00554                 return status;
00555             }
00556 
00557             //  Restore values.
00558             poly2->verts[tri2->v1]->coord[0] = ty1[0];
00559             poly2->verts[tri2->v1]->coord[1] = ty1[1];
00560             poly2->verts[tri2->v1]->coord[2] = ty1[2];
00561             tri2->a = ta; tri2->b = tb; tri2->c = tc; tri2->d = td;
00562 
00563             poly2->verts[tri2->v0]->coord[0] += ta;
00564             poly2->verts[tri2->v0]->coord[1] += tb;
00565             poly2->verts[tri2->v0]->coord[2] += tc;
00566             //  Compute new normal and linecoeff
00567             newplanecoefficients(poly2, tri2);
00568             linecoeff[0] = tri1->c*tri2->b - tri1->b*tri2->c;
00569             linecoeff[1] = tri1->a*tri2->c - tri1->c*tri2->a;
00570             linecoeff[2] = tri1->b*tri2->a - tri1->a*tri2->b;
00571             dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00572                          linecoeff[1]*linecoeff[1] +
00573                          linecoeff[2]*linecoeff[2]);
00574             linecoeff[0] /= dtemp;
00575             linecoeff[1] /= dtemp;
00576             linecoeff[2] /= dtemp;
00577 
00578             status = tri_tri_intersect(tri1,tri2);
00579             if ( status != CUBIT_SUCCESS )
00580             {
00581                 return status;
00582             }
00583 
00584             //  Restore values.
00585             poly2->verts[tri2->v0]->coord[0] = tx1[0];
00586             poly2->verts[tri2->v0]->coord[1] = tx1[1];
00587             poly2->verts[tri2->v0]->coord[2] = tx1[2];
00588             tri2->a = ta; tri2->b = tb; tri2->c = tc; tri2->d = td;
00589             
00590             continue;
00591         }
00592 
00593         if ( do_imprint == true ) {
00594         //  Don't intersect nearly-coplanar triangles.
00595           if ( fabs(tri1->a*tri2->a + tri1->b*tri2->b +tri1->c*tri2->c) > 0.8 ) continue;
00596        }
00597         double dtemp;
00598         dtemp = sqrt(linecoeff[0]*linecoeff[0] +
00599                      linecoeff[1]*linecoeff[1] +
00600                      linecoeff[2]*linecoeff[2]);
00601         linecoeff[0] /= dtemp;
00602         linecoeff[1] /= dtemp;
00603         linecoeff[2] /= dtemp;
00604 
00605         status = tri_tri_intersect(tri1,tri2);
00606         if ( status != CUBIT_SUCCESS )
00607         {
00608             return status;
00609         }
00610 
00611      } // end of loop over poly2
00612 
00613   } // end of loop over poly1
00614 
00615     //delete entire array.
00616   delete [] boxlist;
00617   
00618   return status;
00619 }
00620 
00621 CubitStatus FBIntersect::tri_tri_intersect(FB_Triangle *tri1,
00622               FB_Triangle *tri2)
00623 {
00624 int ret1, ret2, i;
00625 double xc10[3], xc11[3],xc12[3]; // coords for poly1 triangle
00626 double xc20[3], xc21[3],xc22[3]; // coords for poly2 triangle
00627 double d10, d11, d12, d20, d21, d22; // distance of vertex from plane of other triangle
00628 double tt[4];
00629 int edge_vert_type[4];
00630 CubitStatus status;
00631 
00632   int mydebug = 0;
00633   status = CUBIT_SUCCESS;
00634 
00635   tt[0] = tt[1] = tt[2] = tt[3] = CUBIT_DBL_MAX;
00636 //  Is tri1 entirely on one side of tri2?
00637   for ( i = 0; i < 3; i++ ) {
00638      xc10[i] = poly1->verts[tri1->v0]->coord[i];
00639      xc11[i] = poly1->verts[tri1->v1]->coord[i];
00640      xc12[i] = poly1->verts[tri1->v2]->coord[i];
00641    }
00642 
00643    //  distance of each tri1 vert to plane of tri2
00644    d10 = xc10[0]*tri2->a + xc10[1]*tri2->b + xc10[2]*tri2->c + tri2->d;
00645    d11 = xc11[0]*tri2->a + xc11[1]*tri2->b + xc11[2]*tri2->c + tri2->d;
00646    d12 = xc12[0]*tri2->a + xc12[1]*tri2->b + xc12[2]*tri2->c + tri2->d;
00647    if ( ( (d10 < -EPSILON2) && (d11 < -EPSILON2) && (d12 < -EPSILON2) ) ||
00648         ( (d10 > EPSILON2) && (d11 > EPSILON2) && (d12 > EPSILON2) ) ) 
00649       return CUBIT_SUCCESS;
00650 //  Is tri2 entirely on one side of tri1?
00651   for ( i = 0; i < 3; i++ ) {
00652      xc20[i] = poly2->verts[tri2->v0]->coord[i];
00653      xc21[i] = poly2->verts[tri2->v1]->coord[i];
00654      xc22[i] = poly2->verts[tri2->v2]->coord[i];
00655    }
00656 
00657    //  distance of each tri2 vert to plane of tri1
00658    d20 = xc20[0]*tri1->a + xc20[1]*tri1->b + xc20[2]*tri1->c + tri1->d;
00659    d21 = xc21[0]*tri1->a + xc21[1]*tri1->b + xc21[2]*tri1->c + tri1->d;
00660    d22 = xc22[0]*tri1->a + xc22[1]*tri1->b + xc22[2]*tri1->c + tri1->d;
00661    if ( ( (d20 < -EPSILON2) && (d21 < -EPSILON2) && (d22 < -EPSILON2) ) ||
00662         ( (d20 > EPSILON2) && (d21 > EPSILON2) && (d22 > EPSILON2) ) )
00663      return CUBIT_SUCCESS;
00664 //  Get a point on the line of intersection to serve as a reference point.
00665    double ta, tb, ts1, ts2, tdot11;
00666    ts1 = -tri1->d; ts2 = -tri2->d;
00667 
00668    tdot11 = tri1->a*tri2->a + tri1->b*tri2->b + tri1->c*tri2->c;
00669    ta = (ts2*tdot11 - ts1)/(tdot11*tdot11 - 1);
00670    tb = (ts1*tdot11 - ts2)/(tdot11*tdot11 - 1);
00671    linept[0] = ta*tri1->a + tb*tri2->a;
00672    linept[1] = ta*tri1->b + tb*tri2->b;
00673    linept[2] = ta*tri1->c + tb*tri2->c;
00674   
00675 //  There are several cases for the distances.  
00676 //  ret1 holds the number of intersections of the triangle with the
00677 //  intersection line. 
00678 //  Do tri1.
00679       ret1 = get_intersectionline_parameter_values(d10,d11,d12,
00680                                             xc10,xc11,xc12,
00681                    tt[0],tt[1],
00682                    edge_vert_type[0],edge_vert_type[1]);
00683 //  Do tri2.
00684       ret2 = get_intersectionline_parameter_values(d20,d21,d22,
00685                                             xc20,xc21,xc22,
00686                    tt[2],tt[3],
00687                    edge_vert_type[2],edge_vert_type[3]);
00688 //  If not two intersections for each triangle, no intersection edge exists.
00689     if ( (ret1 == 2) && (ret2 == 2) )
00690     {
00691         status = add_intersection_edges(tri1,tri2,tt,edge_vert_type);
00692         if ( status != CUBIT_SUCCESS )
00693         {
00694             return status;
00695         }
00696     }
00697     if(mydebug){
00698       double min_angle_1, max_angle_1;
00699       double min_angle_2, max_angle_2;
00700       poly1->min_max_angles_in_fb_triangle(tri1, min_angle_1, max_angle_1);
00701       poly2->min_max_angles_in_fb_triangle(tri2, min_angle_2, max_angle_2);
00702     
00703     
00704       if(min_angle_1 < 10 || min_angle_2 < 10){
00705         PRINT_INFO(" Tri 1 min angle = %f\n",min_angle_1);
00706         PRINT_INFO(" Tri 2 min angle = %f\n",min_angle_2);
00707       }
00708     }
00709     
00710     return status;
00711 }
00712 
00713 CubitStatus FBIntersect::add_intersection_edges(FB_Triangle *tri1,
00714                    FB_Triangle *tri2,
00715                    double *tt,
00716                    int *edge_vert_type
00717                    )
00718 {
00719 int v10, v11, v20, v21;
00720 bool ifoundit, exists;
00721 FB_Edge *edge;
00722 double x1pt, y1pt, z1pt, x2pt, y2pt, z2pt;
00723 int edge1_vert0_type, edge1_vert1_type, edge2_vert0_type, edge2_vert1_type;
00724 
00725   ifoundit = false;
00726   edge1_vert0_type = edge1_vert1_type = edge2_vert0_type = edge2_vert1_type = UNKNOWN;
00727  
00728 //  Try to handle the epsilon cases by forcing tt[] values that are almost 
00729 //  equal to be equal.
00730 
00731   if ( fabs(tt[0]-tt[2]) < EPSILON ) tt[2] = tt[0];
00732   if ( fabs(tt[0]-tt[3]) < EPSILON ) tt[3] = tt[0];
00733   if ( fabs(tt[1]-tt[2]) < EPSILON ) tt[2] = tt[1];
00734   if ( fabs(tt[1]-tt[3]) < EPSILON ) tt[3] = tt[1];
00735    
00736 // cases 0 and 9, no overlap 
00737 //  if ( (tt[1] < tt[2]) || (tt[3] < tt[0]) ) return CUBIT_SUCCESS;
00738   if ( tt[1] < tt[2] ) { return CUBIT_SUCCESS; }
00739   
00740   if ( tt[3] < tt[0] ) { return CUBIT_SUCCESS; }
00741 
00742 //  These next four cases are by far the most common forms of overlap,
00743 //  so check for them first.
00744 
00745 //  The 1's (e.g.: 11111111) refer to the portion of the line of intersection
00746 //  that is in tri1; the 2's for tri2.  The case numbers were arbirtarily assigned.
00747 
00748 // case 4
00749 //    11111111
00750 //      2222        
00751   if ( (tt[2] > tt[0]) && (tt[3] < tt[1]) ) {
00752     ifoundit = true;    
00753     get_point_from_parameter(tt[2],&x1pt,&y1pt,&z1pt);
00754     get_point_from_parameter(tt[3],&x2pt,&y2pt,&z2pt);
00755     edge1_vert0_type = INTERIOR_VERT;
00756     edge1_vert1_type = INTERIOR_VERT;
00757     edge2_vert0_type = edge_vert_type[2];
00758     edge2_vert1_type = edge_vert_type[3];
00759     
00760 // case 2
00761 //    11111111
00762 //        22222222     
00763   } else if ( (tt[2] < tt[1]) && (tt[2] > tt[0]) && (tt[3] > tt[1]) ) {
00764     ifoundit = true;
00765     get_point_from_parameter(tt[2],&x1pt,&y1pt,&z1pt);
00766     get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);
00767     edge1_vert0_type = INTERIOR_VERT;
00768     edge1_vert1_type = edge_vert_type[1];
00769     edge2_vert0_type = edge_vert_type[2];
00770     edge2_vert1_type = INTERIOR_VERT;
00771   
00772 // case 7
00773 //        11111111
00774 //    22222222      
00775   } else if ( (tt[0] < tt[3]) && (tt[0] > tt[2]) && (tt[1] > tt[3]) ) {
00776     ifoundit = true;  
00777     get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00778     get_point_from_parameter(tt[3],&x2pt,&y2pt,&z2pt);
00779     edge1_vert0_type = edge_vert_type[0];
00780     edge1_vert1_type = INTERIOR_VERT;
00781     edge2_vert0_type = INTERIOR_VERT;
00782     edge2_vert1_type = edge_vert_type[3];
00783     
00784 // case 12
00785 //      1111
00786 //    22222222      
00787   } else if ( (tt[0] > tt[2]) && (tt[1] < tt[3]) ) {
00788     ifoundit = true;  
00789     get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00790     get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);
00791     edge1_vert0_type = edge_vert_type[0];
00792     edge1_vert1_type = edge_vert_type[1];
00793     edge2_vert0_type = INTERIOR_VERT;
00794     edge2_vert1_type = INTERIOR_VERT;
00795     
00796 // case 1 
00797 //    11111111
00798 //           22222222 
00799   } else if ( fabs(tt[1]-tt[2]) < EPSILON ) {
00800     ifoundit = true;
00801 //    return CUBIT_SUCCESS;
00802     get_point_from_parameter(tt[1],&x1pt,&y1pt,&z1pt);
00803     get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);  
00804     edge1_vert0_type = edge_vert_type[1];
00805     edge1_vert1_type = edge_vert_type[1];
00806     edge2_vert0_type = edge_vert_type[2];
00807     edge2_vert1_type = edge_vert_type[2];    
00808    
00809 // case 3
00810   } else if ( tt[0] < tt[2] ) {
00811 // case 3
00812 //    11111111
00813 //       22222    
00814     if ( tt[1] == tt[3] ) {
00815       ifoundit = true;
00816       get_point_from_parameter(tt[2],&x1pt,&y1pt,&z1pt);
00817       get_point_from_parameter(tt[3],&x2pt,&y2pt,&z2pt);  
00818       edge1_vert0_type = INTERIOR_VERT;
00819       edge1_vert1_type = edge_vert_type[1];
00820       edge2_vert0_type = edge_vert_type[2];
00821       edge2_vert1_type = edge_vert_type[3];    
00822     }
00823 
00824 // case 5, 6, or 11   
00825   } else if ( fabs(tt[0]-tt[2]) < EPSILON ) {
00826 
00827 // case 5
00828 //    11111111
00829 //    22222222  
00830     if ( tt[1] == tt[3] ) { 
00831       ifoundit = true;
00832       get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00833       get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);
00834       edge1_vert0_type = edge_vert_type[0];
00835       edge1_vert1_type = edge_vert_type[1];
00836       edge2_vert0_type = edge_vert_type[2];
00837       edge2_vert1_type = edge_vert_type[3];    
00838 
00839 // case 6
00840 //    11111111
00841 //    2222     
00842     } else if ( tt[1] > tt[3] ) { 
00843       ifoundit = true;
00844       get_point_from_parameter(tt[2],&x1pt,&y1pt,&z1pt);
00845       get_point_from_parameter(tt[3],&x2pt,&y2pt,&z2pt);
00846       edge1_vert0_type = edge_vert_type[0];
00847       edge1_vert1_type = INTERIOR_VERT;
00848       edge2_vert0_type = edge_vert_type[2];
00849       edge2_vert1_type = edge_vert_type[3];    
00850 
00851 // case 11
00852 //    1111
00853 //    22222222      
00854     } else {
00855       ifoundit = true;
00856       get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00857       get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);
00858       edge1_vert0_type = edge_vert_type[0];
00859       edge1_vert1_type = edge_vert_type[1];
00860       edge2_vert0_type = edge_vert_type[2];
00861       edge2_vert1_type = INTERIOR_VERT;    
00862     }
00863 
00864 // case 8 or 10      
00865   } else {
00866 
00867 // case 8
00868 //           11111111
00869 //    22222222  
00870     if ( fabs(tt[0]-tt[3]) < EPSILON ) {
00871       ifoundit = true;
00872 //      return CUBIT_SUCCESS;
00873       get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00874       get_point_from_parameter(tt[0],&x2pt,&y2pt,&z2pt);  
00875       edge1_vert0_type = edge_vert_type[0];
00876       edge1_vert1_type = edge_vert_type[0];
00877       edge2_vert0_type = edge_vert_type[3];
00878       edge2_vert1_type = edge_vert_type[3];    
00879 
00880 // case 10
00881 //        1111
00882 //    22222222    
00883     } else if ( fabs(tt[1]-tt[3]) < EPSILON ) {
00884       ifoundit = true;
00885       get_point_from_parameter(tt[0],&x1pt,&y1pt,&z1pt);
00886       get_point_from_parameter(tt[1],&x2pt,&y2pt,&z2pt);       
00887       edge1_vert0_type = edge_vert_type[0];
00888       edge1_vert1_type = edge_vert_type[1];
00889       edge2_vert0_type = INTERIOR_VERT;
00890       edge2_vert1_type = edge_vert_type[3];    
00891     }
00892     
00893   }
00894 
00895   if ( ifoundit == false ) {
00896     PRINT_ERROR("unaccounted for case in add_intersection_edges: tt[] =  %e %e %e %e\n",
00897       tt[0],tt[1],tt[2],tt[3]);
00898     return CUBIT_FAILURE;
00899   } else {
00900     tri1->dudded = true;
00901     v10 = poly1->addavertex(x1pt,y1pt,z1pt);
00902     v11 = poly1->addavertex(x2pt,y2pt,z2pt);
00903     if ( v10 != v11 ) {
00904       exists = poly1->edge_exists_in_tri(*tri1,v10,v11);
00905       if ( exists == false ) {
00906         if ( edge1_vert0_type == INTERIOR_VERT )
00907           edge1_vert0_type = determine_edge_vert_type(edge_vert_type[0],
00908                                                   edge_vert_type[1]);
00909         if ( edge1_vert1_type == INTERIOR_VERT )
00910           edge1_vert1_type = determine_edge_vert_type(edge_vert_type[0],
00911                                                   edge_vert_type[1]);
00912         edge = new FB_Edge(v10,v11,edge1_vert0_type,edge1_vert1_type,true);
00913         tri1->edge_list.push_back(edge);
00914         if ( poly1->edge_exists(v10,v11) == false )
00915           poly1->intersection_edges.push_back(edge);
00916       }
00917     }  else {
00918       int edge_type = UNKNOWN;
00919       int vtype1 = UNKNOWN_VERT, vtype2 = UNKNOWN_VERT;
00920       int v_other1 = UNKNOWN_VERT, v_other2 = UNKNOWN_VERT;
00921 //       edge_type = UNKNOWN;
00922       if ( edge1_vert0_type == EDGE_2 ) {
00923         v_other1 = tri1->v0;
00924         v_other2 = tri1->v2;
00925         edge_type = EDGE_2;
00926         vtype1 = VERTEX_0; 
00927         vtype2 = VERTEX_2;
00928       } else if ( edge1_vert0_type == EDGE_1 ) {
00929         v_other1 = tri1->v1;
00930         v_other2 = tri1->v2;
00931         edge_type = EDGE_1;
00932         vtype1 = VERTEX_1; 
00933         vtype2 = VERTEX_2;
00934       
00935       } else if ( edge1_vert0_type == EDGE_0 ) {
00936         v_other1 = tri1->v0;
00937         v_other2 = tri1->v1;
00938         edge_type = EDGE_0;
00939         vtype1 = VERTEX_0; 
00940         vtype2 = VERTEX_1;
00941       
00942       } 
00943       if ( edge_type != UNKNOWN ) {
00944         exists = poly1->edge_exists_in_tri(*tri1,v_other1,v10);
00945         if ( exists == false ) {    
00946           edge = new FB_Edge(v_other1,v10,vtype1,edge_type,false);
00947           tri1->edge_list.push_back(edge);
00948           if ( poly1->edge_exists(v_other1,v10) == false )
00949             poly1->intersection_edges.push_back(edge);
00950         }
00951         exists = poly1->edge_exists_in_tri(*tri1,v10,v_other2);
00952         if ( exists == false ) {    
00953           edge = new FB_Edge(v10,v_other2,edge_type,vtype2,false);
00954           tri1->edge_list.push_back(edge);
00955           if ( poly1->edge_exists(v10,v_other2) == false )
00956             poly1->intersection_edges.push_back(edge);
00957         } 
00958       }  
00959     }
00960     
00961     tri2->dudded = true;
00962     v20 = poly2->addavertex(x1pt,y1pt,z1pt);
00963     v21 = poly2->addavertex(x2pt,y2pt,z2pt);
00964     if ( v20 != v21) {
00965       exists = poly2->edge_exists_in_tri(*tri2,v20,v21);
00966       if ( exists == false ) {
00967         if ( edge2_vert0_type == INTERIOR_VERT )
00968           edge2_vert0_type = determine_edge_vert_type(edge_vert_type[2],
00969                                                   edge_vert_type[3]);
00970     if ( edge2_vert1_type == INTERIOR_VERT )
00971           edge2_vert1_type = determine_edge_vert_type(edge_vert_type[2],
00972                                                   edge_vert_type[3]);
00973         edge = new FB_Edge(v20,v21,edge2_vert0_type,edge2_vert1_type,true);
00974         tri2->edge_list.push_back(edge);
00975         if ( poly2->edge_exists(v20,v21) == false )
00976           poly2->intersection_edges.push_back(edge);
00977       }
00978     }  else {
00979       int edge_type = UNKNOWN;
00980       int vtype1 = UNKNOWN_VERT, vtype2 = UNKNOWN_VERT;
00981       int v_other1 = UNKNOWN_VERT, v_other2 = UNKNOWN_VERT;
00982 //       edge_type = UNKNOWN;
00983       if ( edge2_vert0_type == EDGE_2 ) {
00984         v_other1 = tri2->v0;
00985         v_other2 = tri2->v2;
00986         edge_type = EDGE_2;
00987         vtype1 = VERTEX_0; 
00988         vtype2 = VERTEX_2;
00989       } else if ( edge2_vert0_type == EDGE_1 ) {
00990         v_other1 = tri2->v1;
00991         v_other2 = tri2->v2;
00992         edge_type = EDGE_1;
00993         vtype1 = VERTEX_1; 
00994         vtype2 = VERTEX_2;
00995       
00996       } else if ( edge2_vert0_type == EDGE_0 ) {
00997         v_other1 = tri2->v0;
00998         v_other2 = tri2->v1;
00999         edge_type = EDGE_0;
01000         vtype1 = VERTEX_0; 
01001         vtype2 = VERTEX_1;
01002       
01003       } 
01004       if ( edge_type != UNKNOWN ) {
01005         exists = poly2->edge_exists_in_tri(*tri2,v_other1,v20);
01006         if ( exists == false ) {    
01007           edge = new FB_Edge(v_other1,v20,vtype1,edge_type,false);
01008           tri2->edge_list.push_back(edge);
01009           if ( poly2->edge_exists(v_other1,v20) == false )
01010             poly2->intersection_edges.push_back(edge);
01011         }
01012         exists = poly2->edge_exists_in_tri(*tri2,v20,v_other2);
01013         if ( exists == false ) {    
01014           edge = new FB_Edge(v20,v_other2,edge_type,vtype2,false);
01015           tri2->edge_list.push_back(edge);
01016           if ( poly2->edge_exists(v20,v_other2) == false )
01017             poly2->intersection_edges.push_back(edge);
01018         }        
01019       }
01020     }
01021   }
01022 
01023   return CUBIT_SUCCESS;
01024   
01025 }
01026 
01027 void FBIntersect::get_point_from_parameter(double parameter,
01028                  double *x, double *y, double *z)
01029 {
01030   *x = linept[0] + linecoeff[0]*parameter; 
01031   *y = linept[1] + linecoeff[1]*parameter; 
01032   *z = linept[2] + linecoeff[2]*parameter; 
01033 
01034 }
01035 
01036 double FBIntersect::get_distance_parameter(double *xc0,
01037                  double *xc1,
01038                  double d0, double d1)
01039 {
01040 double v1dot, v2dot;
01041 
01042 //  dot the coordinates onto the line.
01043 //  Assume that it has been checked already that d0 != d1.
01044 
01045   v1dot = linecoeff[0]*(xc0[0] - linept[0]) +
01046           linecoeff[1]*(xc0[1] - linept[1]) +
01047           linecoeff[2]*(xc0[2] - linept[2]);
01048   v2dot = linecoeff[0]*(xc1[0] - linept[0]) +
01049           linecoeff[1]*(xc1[1] - linept[1]) +
01050           linecoeff[2]*(xc1[2] - linept[2]);
01051 
01052   return v1dot + (v2dot - v1dot)*d0/(d0-d1);
01053 }
01054 
01055 
01056 double FBIntersect::get_distance_parameter_single(double *xc)
01057 {
01058 double coeff;
01059 int coord;
01060 
01061   if ( (fabs(linecoeff[0]) >= fabs(linecoeff[1])) && 
01062        (fabs(linecoeff[0]) >= fabs(linecoeff[2])) ){
01063     coord = 0;
01064   } else if ( fabs(linecoeff[1]) >= fabs(linecoeff[2]) ) {
01065     coord = 1;
01066   }  else {
01067     coord = 2;
01068   }
01069   coeff = linecoeff[coord];
01070   return (xc[coord] - linept[coord])/coeff; 
01071 }
01072 
01073 int FBIntersect::get_intersectionline_parameter_values(
01074                     double d0, 
01075                     double d1, 
01076                     double d2,
01077                     double *pt0,
01078                     double *pt1,
01079                     double *pt2,
01080                     double& t0,
01081                     double& t1,
01082                     int& vert_type_0,
01083                     int& vert_type_1)
01084 {
01085 int ret = 0;
01086 //  ret holds the nnumber of intersections of the line with the triangle.
01087 
01088   if ( fabs(d0) < EPSILON ) {
01089     if ( fabs(d1) < EPSILON ) {
01090       if ( fabs(d2) < EPSILON ) {
01091       // coplanar 
01092         ret = 0;
01093       } else {
01094       //  d0 = d1 = 0
01095         t0 = get_distance_parameter_single(pt0);
01096         t1 = get_distance_parameter_single(pt1);
01097    vert_type_0 = VERTEX_0;
01098    vert_type_1 = VERTEX_1;
01099         ret = 2;
01100       }
01101     } else if ( fabs(d2) < EPSILON ) {
01102     //  d0 = d2 = 0
01103       t0 = get_distance_parameter_single(pt0);
01104       t1 = get_distance_parameter_single(pt2);      
01105       vert_type_0 = VERTEX_0;
01106       vert_type_1 = VERTEX_2;     
01107       ret = 2;
01108     } else if ( d1*d2 < 0.0 ) {
01109     //  d0 = 0 and edge 12 crosses
01110       t0 = get_distance_parameter_single(pt0);
01111       t1 = get_distance_parameter(pt1,pt2,d1,d2);
01112       vert_type_0 = VERTEX_0;
01113       vert_type_1 = EDGE_1;     
01114       ret = 2;
01115     } else {
01116     //  d0 = 0 and no edges cross
01117       ret = 0;
01118     }
01119   } else if ( fabs(d1) < EPSILON ) {
01120     if ( fabs(d2) < EPSILON ) {
01121     //  d1 = d2 = 0
01122       t0 = get_distance_parameter_single(pt1);
01123       t1 = get_distance_parameter_single(pt2); 
01124       vert_type_0 = VERTEX_1;
01125       vert_type_1 = VERTEX_2;     
01126       ret = 2;   
01127     } else if ( d0*d2 < 0.0 ) {
01128     //  d1 = 0 and edge 20 crosses
01129       t0 = get_distance_parameter_single(pt1);    
01130       t1 = get_distance_parameter(pt0,pt2,d0,d2);
01131       vert_type_0 = VERTEX_1;
01132       vert_type_1 = EDGE_2;     
01133       ret = 2;
01134     } else {
01135     //  d1 = 0 and no edges cross
01136       ret = 0;
01137     }
01138   } else if ( fabs(d2) < EPSILON ) {
01139     if ( d0*d1 < 0.0 ) {
01140     //  d2 = 0 and edge 01 crosses
01141       t0 = get_distance_parameter_single(pt2);    
01142       t1 = get_distance_parameter(pt0,pt1,d0,d1);
01143       vert_type_0 = VERTEX_2;
01144       vert_type_1 = EDGE_0;     
01145       ret = 2;    
01146     } else {
01147     //  d2 = 0 and no edges cross
01148       ret = 0;
01149     }
01150   } else if ( d0*d1 < 0.0 ) {
01151     if ( d0*d2 < 0.0 ) {
01152     //  edges 01 and 02 cross
01153       t0 = get_distance_parameter(pt0,pt1,d0,d1);
01154       t1 = get_distance_parameter(pt0,pt2,d0,d2);
01155       vert_type_0 = EDGE_0;
01156       vert_type_1 = EDGE_2;     
01157       ret = 2;
01158     } else {
01159     //  edges 01 and 12 cross
01160       t0 = get_distance_parameter(pt0,pt1,d0,d1);
01161       t1 = get_distance_parameter(pt1,pt2,d1,d2);
01162       vert_type_0 = EDGE_0;
01163       vert_type_1 = EDGE_1;     
01164       ret = 2;
01165     }
01166   } else {
01167   //  edges 02 and 12 cross
01168       t0 = get_distance_parameter(pt0,pt2,d0,d2);
01169       t1 = get_distance_parameter(pt1,pt2,d1,d2);
01170       vert_type_0 = EDGE_2;
01171       vert_type_1 = EDGE_1;     
01172       ret = 2;
01173   }
01174   
01175   if ( ret == 2 ) {  //  Sort them if necessary.
01176     if ( t0 > t1 ) {
01177       int itemp;
01178       double dtemp;
01179       dtemp = t0;
01180       t0 = t1;
01181       t1 = dtemp;
01182       itemp = vert_type_0;
01183       vert_type_0 = vert_type_1;
01184       vert_type_1 = itemp;
01185     }
01186   
01187   } 
01188   
01189   return ret;
01190 }
01191 
01192 void FBIntersect::set_classify_flag(bool value)
01193 {
01194   do_classify = value;
01195 }
01196 
01197 CubitStatus FBIntersect::get_persistent_entity_info(bool *surfs_in,
01198                                                     bool *curves_in,
01199                                                     bool *surfs_out,
01200                                                     bool *curves_out,
01201                                                     const CubitFacetboolOp op,
01202                                                     const int whichparent
01203                                        )
01204 {
01205 std::vector<int> *group, *groupcharacterization;
01206 std::vector<int>::iterator it, ig;
01207 int booltest1, booltest2, booltest1a, booltest2a, booltest_in, booltest_out;
01208 FBPolyhedron *poly;
01209 
01210   if ( nothing_intersected == true ) return CUBIT_SUCCESS;
01211   if ( !poly1 || !poly2 ) {
01212     PRINT_ERROR("Error:  Objects for Booleans must first be created.\n");
01213     return CUBIT_FAILURE;
01214   }
01215   if ( !classify1 || !classify2 ) {
01216     PRINT_ERROR("Error:  Objects for Booleans must first be classified.\n");
01217     return CUBIT_FAILURE;
01218   } 
01219   if ( (whichparent < 1) || (whichparent > 2) ) {
01220     PRINT_ERROR("Error:  Requested nonexistent object.\n");
01221     return CUBIT_FAILURE;  
01222   }
01223   if ( op == CUBIT_FB_UNION ) {
01224     booltest1 = FB_ORIENTATION_INSIDE + FB_ORIENTATION_SAME;
01225     booltest2 = FB_ORIENTATION_INSIDE;
01226     booltest1a = FB_ORIENTATION_OUTSIDE + FB_ORIENTATION_SAME;
01227     booltest2a = FB_ORIENTATION_OUTSIDE;
01228   } else if ( op == CUBIT_FB_INTERSECTION ) {
01229     booltest1 = FB_ORIENTATION_INSIDE + FB_ORIENTATION_SAME;
01230     booltest2 = FB_ORIENTATION_INSIDE; 
01231     booltest1a = FB_ORIENTATION_OUTSIDE+ FB_ORIENTATION_OPPOSITE;
01232     booltest2a = FB_ORIENTATION_INSIDE;         
01233   } else if ( op == CUBIT_FB_SUBTRACTION ) {
01234     booltest1 = FB_ORIENTATION_OUTSIDE+ FB_ORIENTATION_OPPOSITE;
01235     booltest2 = FB_ORIENTATION_INSIDE;    
01236     booltest1a = FB_ORIENTATION_INSIDE + FB_ORIENTATION_SAME;
01237     booltest2a = FB_ORIENTATION_INSIDE; 
01238   } else {
01239     PRINT_ERROR("Error:  Unrecognized Boolean operation.\n");
01240     return CUBIT_FAILURE;
01241   } 
01242  
01243   if ( whichparent == 1 ) {
01244     booltest_in = booltest1;
01245     booltest_out = booltest1a;
01246     classify1->get_group(&group, &groupcharacterization);
01247     poly = poly1;
01248   } else {
01249     booltest_in = booltest2;
01250     booltest_out = booltest2a;
01251     classify2->get_group(&group, &groupcharacterization);
01252     poly = poly2;
01253   }
01254   unsigned int i;
01255   int isurf, icurve0, icurve1, icurve2;
01256 
01257   it = group->begin();
01258   ig = groupcharacterization->begin(); 
01259 
01260   i = 0;
01261   while ( it != group->end() ) {
01262     if ( ( *(ig + *it) & booltest_in ) != 0 ) {
01263 //      vtx = poly1->tris[i]->v0;
01264       isurf = poly->tris[i]->cubitsurfaceindex;
01265       if ( isurf > 0 ) 
01266         surfs_in[isurf] = true;
01267       icurve0 = poly->tris[i]->cubitedge0index;
01268       if ( icurve0 > 0 ) 
01269         curves_in[icurve0] = true;
01270       icurve1 = poly->tris[i]->cubitedge1index;
01271       if ( icurve1 > 0 ) 
01272         curves_in[icurve1] = true;
01273        icurve2 = poly->tris[i]->cubitedge2index;
01274       if ( icurve2 > 0 ) 
01275         curves_in[icurve2] = true;
01276    }
01277     it++;
01278     i++;
01279   }  
01280   
01281   it = group->begin();
01282   ig = groupcharacterization->begin(); 
01283 
01284   i = 0;
01285   while ( it != group->end() ) {
01286     if ( ( *(ig + *it) & booltest_out ) != 0 ) {
01287       isurf = poly->tris[i]->cubitsurfaceindex;
01288       if ( isurf > 0 ) 
01289         surfs_out[isurf] = true;
01290       icurve0 = poly->tris[i]->cubitedge0index;
01291       if ( icurve0 > 0 ) 
01292         curves_out[icurve0] = true;
01293       icurve1 = poly->tris[i]->cubitedge1index;
01294       if ( icurve1 > 0 ) 
01295         curves_out[icurve1] = true;
01296        icurve2 = poly->tris[i]->cubitedge2index;
01297       if ( icurve2 > 0 ) 
01298         curves_out[icurve2] = true;
01299    }
01300     it++;
01301     i++;
01302   }  
01303 
01304   return CUBIT_SUCCESS;
01305 
01306 }
01307                                          
01308 CubitStatus FBIntersect::update_surfs_and_curves(std::vector<double>& out_coords,
01309                                     std::vector<int>& out_connections,
01310                                     std::vector<int> *out_surf_index,
01311                                     std::vector<int> *out_curve_index, 
01312                                     const int whichone)
01313 {
01314 unsigned int i;
01315 FBPolyhedron *poly;
01316  
01317   if ( whichone == 1 ) poly = poly1;
01318   else poly = poly2;
01319   for ( i = 0; i < poly->verts.size(); i++ ) {
01320     out_coords.push_back(poly->verts[i]->coord[0]);
01321     out_coords.push_back(poly->verts[i]->coord[1]);
01322     out_coords.push_back(poly->verts[i]->coord[2]);
01323   } 
01324   for ( i = 0; i < poly->tris.size(); i++ ) {
01325     if ( poly->tris[i]->dudded == true ) continue;
01326     out_connections.push_back(poly->tris[i]->v0);
01327     out_connections.push_back(poly->tris[i]->v1);
01328     out_connections.push_back(poly->tris[i]->v2);
01329     out_surf_index->push_back(poly->tris[i]->cubitsurfaceindex);
01330     out_curve_index->push_back(poly->tris[i]->cubitedge0index);
01331     out_curve_index->push_back(poly->tris[i]->cubitedge1index);
01332     out_curve_index->push_back(poly->tris[i]->cubitedge2index); 
01333   }
01334   return CUBIT_SUCCESS;
01335 }
01336 
01337 CubitStatus FBIntersect::gather_by_boolean(std::vector<double>& out_coords,
01338                                     std::vector<int>& out_connections,
01339                                     std::vector<int> *out_surf_index,
01340                                     std::vector<int> *out_curve_index,
01341                                     std::vector<bool> *is_body_1,
01342                                     const CubitFacetboolOp op)
01343 {
01344     CubitStatus status;
01345     std::vector<int> *group1,
01346                      *group2,
01347                      *groupcharacterization1, 
01348                      *groupcharacterization2;
01349     std::vector<int>::iterator it, ig;
01350     int booltest1, booltest2;
01351 
01352     if ( nothing_intersected == true && op !=CUBIT_FB_UNION)
01353     {
01354         return CUBIT_SUCCESS;
01355     }
01356     if ( !poly1 || !poly2 )
01357     {
01358         PRINT_ERROR("Error:  Objects for Booleans must first be created.\n");
01359         return CUBIT_FAILURE;
01360     }
01361     if ( !classify1 || !classify2 )
01362     {
01363         PRINT_ERROR("Error:  Objects for Booleans must first be classified.\n");
01364         return CUBIT_FAILURE;
01365     } 
01366   
01367     if ( op == CUBIT_FB_UNION )
01368     {
01369         booltest1 = FB_ORIENTATION_OUTSIDE + FB_ORIENTATION_SAME;
01370         booltest2 = FB_ORIENTATION_OUTSIDE;
01371     }
01372     else if ( op == CUBIT_FB_INTERSECTION )
01373     {
01374         booltest1 = FB_ORIENTATION_INSIDE + FB_ORIENTATION_SAME;
01375         booltest2 = FB_ORIENTATION_INSIDE;  
01376     }
01377     else if ( op == CUBIT_FB_SUBTRACTION )
01378     {
01379         booltest1 = FB_ORIENTATION_OUTSIDE+ FB_ORIENTATION_OPPOSITE;
01380         booltest2 = FB_ORIENTATION_INSIDE;    
01381     }
01382     else
01383     {
01384         PRINT_ERROR("Error:  Unrecognized Boolean operation.\n");
01385         return CUBIT_FAILURE;
01386     } 
01387 
01388     classify1->get_group(&group1, &groupcharacterization1);
01389     classify2->get_group(&group2, &groupcharacterization2);
01390 
01391     IntegerHash *hashobj = new IntegerHash(NUMHASHBINS,20);
01392 
01393     it = group1->begin();
01394     ig = groupcharacterization1->begin(); 
01395     unsigned int i;
01396     int vtx, vertnum1, vertnum2, vertnum3, verts_sofar;
01397 
01398     i = 0;
01399     verts_sofar = 0;
01400     while ( it != group1->end() )
01401     {
01402         if ( ( *(ig + *it) & booltest1 ) != 0 )
01403         {
01404             vtx = poly1->tris[i]->v0;
01405             vertnum1 = get_vertex(poly1, vtx, hashobj, out_coords, verts_sofar);
01406 
01407             vtx = poly1->tris[i]->v1;
01408             vertnum2 = get_vertex(poly1, vtx, hashobj, out_coords, verts_sofar);
01409 
01410             vtx = poly1->tris[i]->v2;
01411             vertnum3 = get_vertex(poly1, vtx, hashobj, out_coords, verts_sofar);
01412 
01413             if ( store_connectivity( out_connections,
01414                                      vertnum1,
01415                                      vertnum2,
01416                                      vertnum3 ) != CUBIT_SUCCESS )
01417             {
01418                 return CUBIT_FAILURE;
01419             }
01420             if ( out_surf_index )
01421             {
01422                 out_surf_index->push_back(poly1->tris[i]->cubitsurfaceindex);
01423             }
01424             if ( out_curve_index )
01425             {
01426                 out_curve_index->push_back(poly1->tris[i]->cubitedge0index);
01427                 out_curve_index->push_back(poly1->tris[i]->cubitedge1index);
01428                 out_curve_index->push_back(poly1->tris[i]->cubitedge2index);
01429             }
01430             if ( is_body_1 ) is_body_1->push_back(true);     
01431         }
01432         it++;
01433         i++;
01434     }
01435     it = group2->begin();
01436     ig = groupcharacterization2->begin();
01437     i = 0;
01438     while ( it != group2->end() )
01439     {
01440         if ( ( *(ig + *it) & booltest2 ) != 0 )
01441         {
01442             vtx = poly2->tris[i]->v0;
01443             vertnum1 = get_vertex(poly2, vtx, hashobj, out_coords, verts_sofar);
01444 
01445             vtx = poly2->tris[i]->v1;
01446             vertnum2 = get_vertex(poly2, vtx, hashobj, out_coords, verts_sofar);
01447 
01448             vtx = poly2->tris[i]->v2;
01449             vertnum3 = get_vertex(poly2, vtx, hashobj, out_coords, verts_sofar);
01450 
01451             if ( op == CUBIT_FB_SUBTRACTION ) //  reverse the winding      
01452             {
01453                 if ( store_connectivity( out_connections,
01454                                          vertnum1,
01455                                          vertnum3,
01456                                          vertnum2 ) != CUBIT_SUCCESS )
01457                 {
01458                     return CUBIT_FAILURE;
01459                 }
01460                 if ( out_curve_index )
01461                 {
01462                     out_curve_index->push_back(poly2->tris[i]->cubitedge2index);
01463                     out_curve_index->push_back(poly2->tris[i]->cubitedge1index);
01464                     out_curve_index->push_back(poly2->tris[i]->cubitedge0index);
01465                 }
01466             }
01467             else
01468             {
01469                 if ( store_connectivity( out_connections,
01470                                          vertnum1,
01471                                          vertnum2,
01472                                          vertnum3 ) != CUBIT_SUCCESS )
01473                 {
01474                     return CUBIT_FAILURE;
01475                 }
01476                 if ( out_curve_index )
01477                 {
01478                     out_curve_index->push_back(poly2->tris[i]->cubitedge0index);
01479                     out_curve_index->push_back(poly2->tris[i]->cubitedge1index);
01480                     out_curve_index->push_back(poly2->tris[i]->cubitedge2index); 
01481                 }
01482             }
01483             if ( out_surf_index ) out_surf_index->push_back(poly2->tris[i]->cubitsurfaceindex);
01484             if ( is_body_1 ) is_body_1->push_back(false);
01485         }
01486         it++;
01487         i++;
01488     } 
01489 //  GfxDebug::clear();
01490 //  poly1->debug_draw_boundary_edges(CUBIT_MAGENTA_INDEX);
01491 //  GfxDebug::mouse_xforms();
01492 //  GfxDebug::clear();
01493 //  poly2->debug_draw_boundary_edges(CUBIT_GREEN_INDEX);
01494     delete hashobj;   
01495  
01496     status = CUBIT_SUCCESS;  
01497     return status;
01498 }
01499 
01500 CubitStatus FBIntersect::store_connectivity
01501 (
01502      std::vector<int>& out_connections,
01503      int vertnum1,
01504      int vertnum2,
01505      int vertnum3
01506 )
01507 {
01508     if ( vertnum1 == vertnum2 ||
01509          vertnum2 == vertnum3 ||
01510          vertnum1 == vertnum3 )
01511     {
01512         PRINT_ERROR( "Cannot continue without generating a degenerate facet.\n" );
01513         return CUBIT_FAILURE;
01514     }
01515     out_connections.push_back(vertnum1);   
01516     out_connections.push_back(vertnum2);   
01517     out_connections.push_back(vertnum3); 
01518     return CUBIT_SUCCESS;
01519 }
01520 
01521 int FBIntersect::get_vertex(FBPolyhedron *poly, int vtx, 
01522                             IntegerHash *hashobj,
01523                             std::vector<double>& out_coords,
01524                             int &num_sofar)
01525 {
01526 double xx, yy, zz, xval, yval, zval;
01527 int i, hashvalue, *hasharrayptr, hasharraysize, hptr, ifoundit;
01528 
01529   xx = poly->verts[vtx]->coord[0];
01530   yy = poly->verts[vtx]->coord[1];
01531   zz = poly->verts[vtx]->coord[2];
01532   hashvalue = makeahashvaluefrom_coord(xx,yy,zz);
01533   hasharrayptr = hashobj->getHashBin(hashvalue,&hasharraysize);
01534   ifoundit = -1;
01535   for ( i = 0; i < hasharraysize; i++ ) {
01536     hptr = hasharrayptr[i];
01537     xval = out_coords[3*hptr];
01538     yval = out_coords[3*hptr+1];
01539     zval = out_coords[3*hptr+2];
01540     if ( ( fabs(xval-xx) < EPSILON ) && 
01541     ( fabs(yval-yy) < EPSILON ) &&
01542     ( fabs(zval-zz) < EPSILON ) ) {
01543       ifoundit = hasharrayptr[i];
01544       break;
01545     }
01546   }
01547   if ( ifoundit == -1 ) {
01548     ifoundit = num_sofar;
01549     hashobj->addtoHashList(hashvalue,num_sofar);
01550     out_coords.push_back(xx);
01551     out_coords.push_back(yy);
01552     out_coords.push_back(zz);    
01553     num_sofar++;
01554   }
01555     
01556   return ifoundit;
01557 }
01558                                    
01559 int FBIntersect::makeahashvaluefrom_coord(double x, double y, double z)
01560 {
01561 double mantissasum;
01562 
01563       if ( fabs(x) < 1.e-3 ) x = 0.0;
01564       if ( fabs(y) < 1.e-3 ) y = 0.0;
01565       if ( fabs(z) < 1.e-3 ) z = 0.0;
01566       mantissasum = (int)(10000.0*fabs(x) + 0.5) + 
01567                     (int)(10000.0*fabs(y) + 0.5) + 
01568           (int)(10000.0*fabs(z) + 0.5);
01569       
01570       return (int)(mantissasum) % NUMHASHBINS;
01571 }
01572                                     
01573                                       
01574 void FBIntersect::set_body1_planar()
01575 {
01576   body1_is_plane = true;
01577 }
01578                               
01579 void FBIntersect::set_body2_planar()
01580 {
01581   body2_is_plane = true;
01582 }
01583 
01584 void FBIntersect::set_imprint()
01585 {
01586   do_imprint = true;
01587 }
01588 
01589 void FBIntersect::newplanecoefficients(FBPolyhedron *poly, FB_Triangle *tri)
01590 {
01591 FB_Coord *mycoord;
01592 double x1, x2, x3, y1, y2, y3, z1, z2, z3, e1x, e1y, e1z, e2x, e2y, e2z;
01593 double a, b, c, d, dtemp;
01594 
01595      mycoord = poly->verts[tri->v0];
01596      x1 = mycoord->coord[0];
01597      y1 = mycoord->coord[1];
01598      z1 = mycoord->coord[2];
01599      mycoord = poly->verts[tri->v1];
01600      x2 = mycoord->coord[0];
01601      y2 = mycoord->coord[1];
01602      z2 = mycoord->coord[2];
01603      mycoord = poly->verts[tri->v2];
01604      x3 = mycoord->coord[0];
01605      y3 = mycoord->coord[1];
01606      z3 = mycoord->coord[2];
01607      e1x = x1 - x2; e1y = y1 - y2; e1z = z1 - z2;
01608      e2x = x3 - x2; e2y = y3 - y2; e2z = z3 - z2;
01609      a = e1z*e2y - e2z*e1y;
01610      b = e1x*e2z - e2x*e1z;
01611      c = e1y*e2x - e2y*e1x;
01612      dtemp = sqrt(a*a + b*b + c*c);
01613      if ( dtemp > EPSILON2 ) {
01614        a /= dtemp;
01615        b /= dtemp;
01616        c /= dtemp;
01617      } else {
01618        PRINT_WARNING("small-area triangle\n");
01619      }
01620      d = -(a*x1 + b*y1 + c*z1);
01621      tri->a = a; tri->b = b; tri->c = c; tri->d = d;
01622           
01623 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines