cgma
FBImprint.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 <math.h>
00029 #include "FBImprint.hpp"
00030 #include "FBPolyhedron.hpp"
00031 #include "FBRetriangulate.hpp"
00032 #include "GeometryDefines.h"
00033 #include "CubitMessage.hpp"
00034 #include "FBDataUtil.hpp"
00035 #include <stdio.h>
00036 
00037 FBImprint::FBImprint()
00038 {
00039 
00040 }
00041 
00042 FBImprint::~FBImprint()
00043 {
00044 
00045 }  
00046 
00047 CubitStatus FBImprint::imprint_body_curve(const std::vector<double>& Bodycoords,
00048                            const std::vector<int>& Bodyconnections,
00049                            const std::vector<FB_Coord*>& FB_imprint_edge_coords,
00050                            const std::vector<FB_Edge*>& FB_imprint_edges,
00051                            const std::vector<FSBoundingBox*>& FB_imprint_edge_bboxes,
00052                            std::vector<int>* indices)
00053 {
00054 CubitStatus status;
00055 bool new_edge_created;
00056 
00057   status = CUBIT_SUCCESS;   
00058 
00059   if ( Bodycoords.size()%3 != 0 ) {
00060     PRINT_ERROR("Bad coordinates for first part fed to FBImprint.\n");
00061     return CUBIT_FAILURE;
00062   }
00063   if ( Bodyconnections.size()%3 != 0 ) {
00064     PRINT_ERROR("Bad connection list for first part fed to FBImprint.\n");
00065     return CUBIT_FAILURE;
00066   }
00067  
00068   f_c_indices = indices;
00069   poly = new FBPolyhedron;
00070   status = poly->makepoly(Bodycoords,Bodyconnections,f_c_indices);
00071   
00072   status = edges_tri_intersect(FB_imprint_edge_coords,FB_imprint_edges,
00073                                FB_imprint_edge_bboxes,new_edge_created);
00074 
00075   if ( new_edge_created == true ) {
00076     std::vector<int> newFacets;
00077     status = poly->retriangulate(newFacets);
00078 /*    
00079     FILE *out;
00080     out = fopen("1Qaz.fct","w");
00081     int numtris, numverts;
00082     numverts = poly->verts.size();
00083 
00084     int ii;  
00085 
00086     numtris = 0;
00087     for ( ii = 0; ii < poly->tris.size(); ii++ ) 
00088       if ( poly->tris[ii]->dudded == false ) numtris++;
00089     fprintf(out,"%d %d\n",numverts,numtris);
00090     for ( ii = 0; ii < poly->verts.size(); ii++ ) 
00091       fprintf(out,"%d %le %le %le\n",ii+1,poly->verts[ii]->coord[0],
00092            poly->verts[ii]->coord[1],poly->verts[ii]->coord[2]);
00093     for ( ii = 0; ii < poly->tris.size(); ii++ ) {
00094       if ( poly->tris[ii]->dudded == false )
00095         fprintf(out,"%d %d %d %d\n",ii+1,1+poly->tris[ii]->v0,  
00096            1+poly->tris[ii]->v1,1+poly->tris[ii]->v2);
00097     }   
00098     fclose(out); 
00099 */        
00100   }  
00101   
00102   return status;
00103 
00104 }
00105 
00106 CubitStatus FBImprint::edges_tri_intersect(const std::vector<FB_Coord*>& FB_imprint_edge_coords,
00107                            const std::vector<FB_Edge*>& FB_imprint_edges,
00108                            const std::vector<FSBoundingBox*>& FB_imprint_edge_bboxes,
00109                            bool &new_edge_created)
00110 {
00111   CubitStatus status;
00112   unsigned int i, j;
00113   int numboxesfound, *boxlist;
00114   FSBoundingBox* edgebox;
00115   double edge_dir[3] = {0.0}, edge_0[3] = {0.0}, edge_1[3] = {0.0}, edge_length = 0.0;
00116   bool big_angle;
00117 
00118   boxlist = new int[poly->tris.size()];
00119 
00120   status = CUBIT_SUCCESS;
00121   new_edge_created = false;
00122   for ( i = 0; i < FB_imprint_edge_bboxes.size(); i++ ) {
00123      edgebox = FB_imprint_edge_bboxes[i];
00124      if ( (edgebox->xmax < poly->polyxmin) || 
00125           (edgebox->xmin > poly->polyxmax) ||
00126           (edgebox->ymax < poly->polyymin) || 
00127           (edgebox->ymin > poly->polyymax) ||
00128           (edgebox->zmax < poly->polyzmin) || 
00129           (edgebox->zmin > poly->polyzmax) ) continue;
00130      poly->kdtree->box_kdtree_intersect(*edgebox,&numboxesfound,boxlist);
00131      if ( numboxesfound > 0 ) {  //  Get a unit vector along the edge.
00132        edge_0[0] = FB_imprint_edge_coords[FB_imprint_edges[i]->v0]->coord[0];
00133        edge_1[0] = FB_imprint_edge_coords[FB_imprint_edges[i]->v1]->coord[0];
00134        edge_dir[0] = edge_1[0] - edge_0[0];
00135        edge_0[1] = FB_imprint_edge_coords[FB_imprint_edges[i]->v0]->coord[1];
00136        edge_1[1] = FB_imprint_edge_coords[FB_imprint_edges[i]->v1]->coord[1];
00137        edge_dir[1] = edge_1[1] - edge_0[1];
00138        edge_0[2] = FB_imprint_edge_coords[FB_imprint_edges[i]->v0]->coord[2];
00139        edge_1[2] = FB_imprint_edge_coords[FB_imprint_edges[i]->v1]->coord[2];
00140        edge_dir[2] = edge_1[2] - edge_0[2];
00141        edge_length = sqrt(edge_dir[0]*edge_dir[0] + edge_dir[1]*edge_dir[1] + edge_dir[2]*edge_dir[2]);
00142        if ( edge_length < GEOMETRY_RESABS ) continue;
00143        edge_dir[0] /= edge_length;
00144        edge_dir[1] /= edge_length;
00145        edge_dir[2] /= edge_length;             
00146      }
00147 
00148      for ( j = 0; j < (unsigned int)numboxesfound; j++ ) { 
00149         FB_Triangle *tri = poly->tris[boxlist[j]];     
00150         if ( (edgebox->xmax < tri->boundingbox.xmin) || 
00151              (edgebox->xmin > tri->boundingbox.xmax) ||
00152              (edgebox->ymax < tri->boundingbox.ymin) || 
00153              (edgebox->ymin > tri->boundingbox.ymax) ||
00154              (edgebox->zmax < tri->boundingbox.zmin) || 
00155              (edgebox->zmin > tri->boundingbox.zmax) ) continue;
00156 //  Try to get a reasonable value for imprint_res.  It should depend on the size of the 
00157 //  triangle and the length of the edge.
00158         double tri_size = sqrt( (tri->boundingbox.xmax-tri->boundingbox.xmin)*
00159                                 (tri->boundingbox.xmax-tri->boundingbox.xmin) +
00160                                 (tri->boundingbox.ymax-tri->boundingbox.ymin)*
00161                                 (tri->boundingbox.ymax-tri->boundingbox.ymin) +
00162                                 (tri->boundingbox.zmax-tri->boundingbox.zmin)*
00163                                 (tri->boundingbox.zmax-tri->boundingbox.zmin) );                                
00164         if ( edge_length < tri_size ) imprint_res = 0.01*edge_length; 
00165         else imprint_res = 0.01*tri_size;          
00166 //  Flag triangles with planes greater than ~15 degree angle wrt the edge.
00167         big_angle = false;
00168         if ( fabs(edge_dir[0]*tri->a +edge_dir[1]*tri->b +edge_dir[2]*tri->c) > 0.25 ) 
00169           big_angle = true;  
00170         status = single_edge_tri_intersect(edge_0,edge_1,new_edge_created,tri,big_angle);
00171      }
00172 
00173   }
00174 
00175   return status;
00176   
00177 }
00178 
00179 //  edge_0 and edge_1 are edge end-points; edge_dir is unit vector from 0 to 1                           
00180 CubitStatus FBImprint::single_edge_tri_intersect(double *edge_0,double *edge_1,
00181                                                  bool &new_edge_created,
00182                                                  FB_Triangle *tri,
00183                                                  bool big_angle) 
00184 {
00185 CubitStatus status;
00186 double distance0, distance1;
00187 
00188   status = CUBIT_SUCCESS;
00189   distance0 = edge_0[0]*tri->a + edge_0[1]*tri->b + edge_0[2]*tri->c + tri->d;                                                                     
00190   distance1 = edge_1[0]*tri->a + edge_1[1]*tri->b + edge_1[2]*tri->c + tri->d;                                                                     
00191 
00192   //  If both end-points are farther away from the plane of the triangle than
00193   //  imprint_res and on the same side, there is no intersection.
00194   if ( ( (distance0 > imprint_res) && (distance1 > imprint_res) ) ||
00195        ( (distance0 < -imprint_res) && (distance1 < -imprint_res) ) ) 
00196     return status;
00197 
00198   //  Check the edge-triangle border closest distances.
00199 
00200   double d0[3], d1[3], s, t, sunclipped, tunclipped;
00201   double closest_dist, tri_pt[3];
00202   bool parallel0, parallel1, parallel2;
00203   int numptsfound = 0;
00204   double edge_intersection_pt[3][2]; 
00205   int edge_vert_type[2];
00206   
00207   //  Test for closest distance from the edge to each of the triangle edges.
00208   //  In order for an edge to be generated in the triangle, the intersection
00209   //  parameter tunclipped has to lie between 0 and 1.  If this condition is met,
00210   //  the next requirement is that the intersection parameter for the test
00211   //  edge, sunclipped, has to be between 0 and 1 and the intersection distance
00212   //  has to be less than imprint_res, or the intersection edge endpoint 
00213   //  distance (distance0 or distance1) has to be less than imprint_res.
00214   
00215   d0[0] = edge_1[0] - edge_0[0];
00216   d0[1] = edge_1[1] - edge_0[1];
00217   d0[2] = edge_1[2] - edge_0[2];
00218   d1[0] = poly->verts[tri->v1]->coord[0] - poly->verts[tri->v0]->coord[0];
00219   d1[1] = poly->verts[tri->v1]->coord[1] - poly->verts[tri->v0]->coord[1];
00220   d1[2] = poly->verts[tri->v1]->coord[2] - poly->verts[tri->v0]->coord[2];
00221   tri_pt[0] = poly->verts[tri->v0]->coord[0];
00222   tri_pt[1] = poly->verts[tri->v0]->coord[1];
00223   tri_pt[2] = poly->verts[tri->v0]->coord[2];
00224   
00225   closest_dist = FBDataUtil::closest_seg_seg_dist(edge_0,d0,tri_pt,d1,&s,&t,
00226               &sunclipped,&tunclipped,&parallel0);
00227               
00228   if ( (tunclipped >= 0.0) && (tunclipped <= 1.0) ) {
00229     if ( (sunclipped >= 0.0) && (sunclipped <= 1.0) && 
00230          (closest_dist <  imprint_res) ) {
00231       edge_intersection_pt[0][numptsfound] = tri_pt[0] + tunclipped*d1[0];
00232       edge_intersection_pt[1][numptsfound] = tri_pt[1] + tunclipped*d1[1];
00233       edge_intersection_pt[2][numptsfound] = tri_pt[2] + tunclipped*d1[2];
00234       if ( tunclipped == 0.0 ) edge_vert_type[numptsfound] = VERTEX_0;
00235       else if ( tunclipped == 1.0 ) edge_vert_type[numptsfound] = VERTEX_1;
00236       else edge_vert_type[numptsfound] = EDGE_0;
00237       numptsfound++;     
00238     }
00239     else if ( sunclipped < 0.0 ) {
00240       if ( fabs(distance0) < imprint_res ) {
00241         edge_intersection_pt[0][numptsfound] = edge_0[0] - distance0*tri->a;
00242         edge_intersection_pt[1][numptsfound] = edge_0[1] - distance0*tri->b;
00243         edge_intersection_pt[2][numptsfound] = edge_0[2] - distance0*tri->c;
00244         edge_vert_type[numptsfound] = INTERIOR_VERT;
00245         numptsfound++;
00246       } 
00247     }
00248     else if ( sunclipped > 1.0 ) {
00249       if ( fabs(distance1) < imprint_res ) {
00250         edge_intersection_pt[0][numptsfound] = edge_1[0] - distance1*tri->a;
00251         edge_intersection_pt[1][numptsfound] = edge_1[1] - distance1*tri->b;
00252         edge_intersection_pt[2][numptsfound] = edge_1[2] - distance1*tri->c;
00253         edge_vert_type[numptsfound] = INTERIOR_VERT;
00254         numptsfound++;      
00255       }    
00256     }
00257   }             
00258               
00259   d1[0] = poly->verts[tri->v2]->coord[0] - poly->verts[tri->v1]->coord[0];
00260   d1[1] = poly->verts[tri->v2]->coord[1] - poly->verts[tri->v1]->coord[1];
00261   d1[2] = poly->verts[tri->v2]->coord[2] - poly->verts[tri->v1]->coord[2];
00262   tri_pt[0] = poly->verts[tri->v1]->coord[0];
00263   tri_pt[1] = poly->verts[tri->v1]->coord[1];
00264   tri_pt[2] = poly->verts[tri->v1]->coord[2];
00265   
00266   closest_dist = FBDataUtil::closest_seg_seg_dist(edge_0,d0,tri_pt,d1,&s,&t,
00267               &sunclipped,&tunclipped,&parallel1);
00268 
00269   if ( (tunclipped >= 0.0) && (tunclipped <= 1.0) ) {
00270     if ( (sunclipped >= 0.0) && (sunclipped <= 1.0) && 
00271          (closest_dist <  imprint_res) ) {
00272       edge_intersection_pt[0][numptsfound] = tri_pt[0] + tunclipped*d1[0];
00273       edge_intersection_pt[1][numptsfound] = tri_pt[1] + tunclipped*d1[1];
00274       edge_intersection_pt[2][numptsfound] = tri_pt[2] + tunclipped*d1[2]; 
00275       if ( tunclipped == 0.0 ) edge_vert_type[numptsfound] = VERTEX_1;
00276       else if ( tunclipped == 1.0 ) edge_vert_type[numptsfound] = VERTEX_2;
00277       else edge_vert_type[numptsfound] = EDGE_1;
00278       numptsfound++;     
00279     }
00280     else if ( sunclipped < 0.0 ) {
00281       if ( fabs(distance0) < imprint_res ) {
00282         edge_intersection_pt[0][numptsfound] = edge_0[0] - distance0*tri->a;
00283         edge_intersection_pt[1][numptsfound] = edge_0[1] - distance0*tri->b;
00284         edge_intersection_pt[2][numptsfound] = edge_0[2] - distance0*tri->c;
00285         edge_vert_type[numptsfound] = INTERIOR_VERT;
00286         numptsfound++;
00287       } 
00288     }
00289     else if ( sunclipped > 1.0 ) {
00290       if ( fabs(distance1) < imprint_res ) {
00291         edge_intersection_pt[0][numptsfound] = edge_1[0] - distance1*tri->a;
00292         edge_intersection_pt[1][numptsfound] = edge_1[1] - distance1*tri->b;
00293         edge_intersection_pt[2][numptsfound] = edge_1[2] - distance1*tri->c;
00294         edge_vert_type[numptsfound] = INTERIOR_VERT;
00295         numptsfound++;      
00296       }    
00297     }
00298   }                          
00299          
00300   if ( numptsfound < 2 ) {         
00301     d1[0] = poly->verts[tri->v0]->coord[0] - poly->verts[tri->v2]->coord[0];
00302     d1[1] = poly->verts[tri->v0]->coord[1] - poly->verts[tri->v2]->coord[1];
00303     d1[2] = poly->verts[tri->v0]->coord[2] - poly->verts[tri->v2]->coord[2];
00304     tri_pt[0] = poly->verts[tri->v2]->coord[0];
00305     tri_pt[1] = poly->verts[tri->v2]->coord[1];
00306     tri_pt[2] = poly->verts[tri->v2]->coord[2];
00307 
00308     closest_dist = FBDataUtil::closest_seg_seg_dist(edge_0,d0,tri_pt,d1,&s,&t,
00309                 &sunclipped,&tunclipped,&parallel2);
00310 
00311     if ( (tunclipped >= 0.0) && (tunclipped <= 1.0) ) {
00312       if ( (sunclipped >= 0.0) && (sunclipped <= 1.0) && 
00313            (closest_dist <  imprint_res) ) {
00314         edge_intersection_pt[0][numptsfound] = tri_pt[0] + tunclipped*d1[0];
00315         edge_intersection_pt[1][numptsfound] = tri_pt[1] + tunclipped*d1[1];
00316         edge_intersection_pt[2][numptsfound] = tri_pt[2] + tunclipped*d1[2]; 
00317         if ( tunclipped == 0.0 ) edge_vert_type[numptsfound] = VERTEX_2;
00318         else if ( tunclipped == 1.0 ) edge_vert_type[numptsfound] = VERTEX_0;
00319         else edge_vert_type[numptsfound] = EDGE_2;
00320         numptsfound++;     
00321       }
00322       else if ( sunclipped < 0.0 ) {
00323         if ( fabs(distance0) < imprint_res ) {
00324           edge_intersection_pt[0][numptsfound] = edge_0[0] - distance0*tri->a;
00325           edge_intersection_pt[1][numptsfound] = edge_0[1] - distance0*tri->b;
00326           edge_intersection_pt[2][numptsfound] = edge_0[2] - distance0*tri->c;
00327           edge_vert_type[numptsfound] = INTERIOR_VERT;
00328           numptsfound++;
00329         } 
00330       }
00331       else if ( sunclipped > 1.0 ) {
00332         if ( fabs(distance1) < imprint_res ) {
00333           edge_intersection_pt[0][numptsfound] = edge_1[0] - distance1*tri->a;
00334           edge_intersection_pt[1][numptsfound] = edge_1[1] - distance1*tri->b;
00335           edge_intersection_pt[2][numptsfound] = edge_1[2] - distance1*tri->c;
00336           edge_vert_type[numptsfound] = INTERIOR_VERT;
00337           numptsfound++;      
00338         }    
00339       }
00340     }
00341   }  
00342   FB_Edge *edge;
00343   int v10, v11;
00344   bool exists;
00345 
00346   if ( numptsfound == 2 ) {
00347     
00348     tri->dudded = true;
00349     v10 = poly->addavertex(edge_intersection_pt[0][0],
00350                            edge_intersection_pt[1][0],
00351                            edge_intersection_pt[2][0]);
00352     v11 = poly->addavertex(edge_intersection_pt[0][1],
00353                            edge_intersection_pt[1][1],
00354                            edge_intersection_pt[2][1]);
00355                     
00356     if ( v10 != v11 ) {
00357       exists = poly->edge_exists_in_tri(*tri,v10,v11);
00358       if ( exists == false ) {
00359         new_edge_created = true;
00360         edge = new FB_Edge(v10,v11,edge_vert_type[0],edge_vert_type[1],true);
00361         tri->edge_list.push_back(edge);
00362         if ( poly->edge_exists(v10,v11) == false )
00363           poly->intersection_edges.push_back(edge);
00364       }
00365     }
00366   } else if ( numptsfound == 1 ) {
00367   //  Is it on an edge?
00368     int edge_type = UNKNOWN;
00369     int vtype1 = UNKNOWN_VERT, vtype2 = UNKNOWN_VERT;
00370     int v_other1 = UNKNOWN_VERT, v_other2 = UNKNOWN_VERT;
00371 //     edge_type = UNKNOWN;                      
00372     if ( edge_vert_type[0] == EDGE_0 ) {
00373         v_other1 = tri->v0;
00374         v_other2 = tri->v1;
00375         edge_type = EDGE_0;
00376         vtype1 = VERTEX_0; 
00377         vtype2 = VERTEX_1;    
00378     } else if ( edge_vert_type[0] == EDGE_1 ) {
00379         v_other1 = tri->v1;
00380         v_other2 = tri->v2;
00381         edge_type = EDGE_1;
00382         vtype1 = VERTEX_1; 
00383         vtype2 = VERTEX_2;
00384     } else if ( edge_vert_type[0] == EDGE_2 ) {
00385         v_other1 = tri->v0;
00386         v_other2 = tri->v2;
00387         edge_type = EDGE_2;
00388         vtype1 = VERTEX_0; 
00389         vtype2 = VERTEX_2;
00390     }
00391     if ( edge_type != UNKNOWN ) {
00392       tri->dudded = true;
00393       v10 = poly->addavertex(edge_intersection_pt[0][0],
00394                              edge_intersection_pt[1][0],
00395                              edge_intersection_pt[2][0]);
00396                              
00397         exists = poly->edge_exists_in_tri(*tri,v_other1,v10);
00398         if ( exists == false ) {    
00399           new_edge_created = true;
00400           edge = new FB_Edge(v_other1,v10,vtype1,edge_type,false);
00401           tri->edge_list.push_back(edge);
00402           if ( poly->edge_exists(v_other1,v10) == false )
00403             poly->intersection_edges.push_back(edge);
00404         }
00405         exists = poly->edge_exists_in_tri(*tri,v10,v_other2);
00406         if ( exists == false ) {    
00407           new_edge_created = true;
00408           edge = new FB_Edge(v10,v_other2,edge_type,vtype2,false);
00409           tri->edge_list.push_back(edge);
00410           if ( poly->edge_exists(v10,v_other2) == false )
00411             poly->intersection_edges.push_back(edge);
00412         } 
00413      
00414     
00415      }
00416   }
00417 
00418   return status;
00419 }                           
00420                                                                     
00421 CubitStatus FBImprint::update_surfs_and_curves(std::vector<double>& out_coords,
00422                                     std::vector<int>& out_connections,
00423                                     std::vector<int> *out_surf_index,
00424                                     std::vector<int> *out_curve_index)
00425 {
00426 unsigned int i;
00427 
00428   for ( i = 0; i < poly->verts.size(); i++ ) {
00429     out_coords.push_back(poly->verts[i]->coord[0]);
00430     out_coords.push_back(poly->verts[i]->coord[1]);
00431     out_coords.push_back(poly->verts[i]->coord[2]);
00432   } 
00433   for ( i = 0; i < poly->tris.size(); i++ ) {
00434     if ( poly->tris[i]->dudded == true ) continue;
00435     out_connections.push_back(poly->tris[i]->v0);
00436     out_connections.push_back(poly->tris[i]->v1);
00437     out_connections.push_back(poly->tris[i]->v2);
00438     out_surf_index->push_back(poly->tris[i]->cubitsurfaceindex);
00439     out_curve_index->push_back(poly->tris[i]->cubitedge0index);
00440     out_curve_index->push_back(poly->tris[i]->cubitedge1index);
00441     out_curve_index->push_back(poly->tris[i]->cubitedge2index);
00442   }
00443   return CUBIT_SUCCESS;
00444 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines