cgma
FBDataUtil.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 "FBDataUtil.hpp"
00028 #include "GeometryDefines.h"
00029 #include "DLIList.hpp"
00030 
00031 //===========================================================================
00032 //Function Name: intersect_plane_with_boundingbox
00033 //Description:   Find and return the number of intersections of a plane with
00034 //               the edges of a bounding box.  There will be 6 or fewer.
00035 //               These will be returned in the intersection_points array.
00036 //Author: mlstate - rewritten from John Fowler's original version.
00037 //Date: 3/2005
00038 //===========================================================================
00039 void FBDataUtil::intersect_plane_with_boundingbox(CubitBox &bbox, 
00040                                      const CubitVector &v0, 
00041                                      const CubitVector &v1,
00042                                      const CubitVector &v2,
00043                                      DLIList<CubitVector> &intersection_points )
00044 {
00045 CubitVector v3 = v0 - v1;
00046 CubitVector v4 = v2 - v1;
00047 CubitVector threepointnormal = v3*v4;
00048 threepointnormal.normalize();
00049 //CubitVector threepointcenter = (v0 + v1 + v2)/3.;
00050 //  Need to get the size and location of the facetted cutting plane.  Do this by 
00051 //  intersecting the 3-point plane with the bounding box edges.  Will get a max of
00052 //  six intersections.  The size of the box for these intersections will determine
00053 //  the cutting plane size; the centroid of the intersections will be the location
00054 //  of the center of the cutting plane.
00055 
00056 double a, b, c, d;  //  plane coefficients
00057   a = threepointnormal.x();
00058   b = threepointnormal.y();
00059   c = threepointnormal.z();
00060   d = -(a*v0.x() + b*v0.y() + c*v0.z());
00061 double xbmin, ybmin, zbmin, xbmax, ybmax, zbmax;
00062 double testpoint;
00063 
00064     CubitVector boxmin = bbox.minimum();
00065     xbmin = boxmin.x(); ybmin = boxmin.y(); zbmin = boxmin.z();
00066     CubitVector boxmax = bbox.maximum();
00067     xbmax = boxmax.x(); ybmax = boxmax.y(); zbmax = boxmax.z();
00068 
00069     if ( fabs(a) < GEOMETRY_RESABS &&
00070          fabs(b) < GEOMETRY_RESABS )
00071     {
00072         // Normal points in Z dir only.
00073         //  a = 0; b = 0
00074         if ( ((zbmin + d/c) < -GEOMETRY_RESABS) && ((zbmax + d/c) > GEOMETRY_RESABS) )
00075         {
00076             add_unique_point( intersection_points, CubitVector(xbmin, ybmin, -d/c) );
00077             add_unique_point( intersection_points, CubitVector(xbmax, ybmin, -d/c) );
00078             add_unique_point( intersection_points, CubitVector(xbmax, ybmax, -d/c) );
00079             add_unique_point( intersection_points, CubitVector(xbmin, ybmax, -d/c) );
00080         }
00081     }
00082     else if ( fabs(a) < GEOMETRY_RESABS &&
00083               fabs(c) < GEOMETRY_RESABS )
00084     {
00085         // Normal points in Y dir only.
00086         if ( ((ybmin + d/b) < -GEOMETRY_RESABS) && ((ybmax + d/b) > GEOMETRY_RESABS) )
00087         {
00088             add_unique_point( intersection_points, CubitVector(xbmin, -d/b, zbmin ) );
00089             add_unique_point( intersection_points, CubitVector(xbmax, -d/b, zbmin ) );
00090             add_unique_point( intersection_points, CubitVector(xbmax, -d/b, zbmax ) );
00091             add_unique_point( intersection_points, CubitVector(xbmin, -d/b, zbmax ) );
00092         }
00093     }
00094     else if ( fabs(b) < GEOMETRY_RESABS &&
00095               fabs(c) < GEOMETRY_RESABS )
00096     {
00097         // Normal points in X dir only.
00098         if ( ((xbmin + d/a) < -GEOMETRY_RESABS) && ((xbmax + d/a) > GEOMETRY_RESABS) )
00099         {
00100             add_unique_point( intersection_points, CubitVector(-d/a, ybmin, zbmin ) );
00101             add_unique_point( intersection_points, CubitVector(-d/a, ybmax, zbmin ) );
00102             add_unique_point( intersection_points, CubitVector(-d/a, ybmax, zbmax ) );
00103             add_unique_point( intersection_points, CubitVector(-d/a, ybmin, zbmax ) );
00104         }
00105     }
00106     else if ( fabs(a) < GEOMETRY_RESABS )
00107     {
00108         // Normal is in YZ plane.
00109         testpoint = -(c*zbmin + d)/b;
00110         if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
00111              (ybmax > (testpoint - GEOMETRY_RESABS)) )
00112         {
00113             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
00114             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
00115         }
00116         testpoint = -(c*zbmax + d)/b;
00117         if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
00118              (ybmax > (testpoint - GEOMETRY_RESABS)) )
00119         {
00120             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
00121             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
00122         }
00123         testpoint = -(b*ybmin + d)/c;
00124         if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
00125              (zbmax > (testpoint - GEOMETRY_RESABS)) )
00126         {
00127             add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
00128             add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
00129         }
00130         testpoint = -(b*ybmax + d)/c;
00131         if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
00132              (zbmax > (testpoint - GEOMETRY_RESABS)) )
00133         {
00134             add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
00135             add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
00136         }
00137     }
00138     else if ( fabs(b) < GEOMETRY_RESABS )
00139     {
00140         // Normal is in XZ plane
00141         testpoint = -(c*zbmin + d)/a;
00142         if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
00143              (xbmax > (testpoint - GEOMETRY_RESABS)) )
00144         {
00145             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
00146             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
00147         }
00148         testpoint = -(c*zbmax + d)/a;
00149         if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
00150              (xbmax > (testpoint - GEOMETRY_RESABS)) )
00151         {
00152             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
00153             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
00154         }
00155         testpoint = -(a*xbmin + d)/c;
00156         if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
00157              (zbmax > (testpoint - GEOMETRY_RESABS)) )
00158         {
00159             add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
00160             add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
00161         }
00162         testpoint = -(a*xbmax + d)/c;
00163         if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
00164              (zbmax > (testpoint - GEOMETRY_RESABS)) )
00165         {
00166             add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
00167             add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
00168         }
00169     }
00170     else if ( fabs(c) < GEOMETRY_RESABS )
00171     {
00172         // Normal is in XY plane
00173         testpoint = -(a*xbmin + d)/b;
00174         if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
00175              (ybmax > (testpoint - GEOMETRY_RESABS)) )
00176         {
00177             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
00178             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
00179         }
00180         testpoint = -(a*xbmax + d)/b;
00181         if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
00182              (ybmax > (testpoint - GEOMETRY_RESABS)) )
00183         {
00184             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
00185             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
00186         }
00187         testpoint = -(b*ybmin + d)/a;
00188         if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
00189              (xbmax > (testpoint - GEOMETRY_RESABS)) )
00190         {
00191             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
00192             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
00193         }
00194         testpoint = -(b*ybmax + d)/a;
00195         if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
00196              (xbmax > (testpoint - GEOMETRY_RESABS)) )
00197         {
00198             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
00199             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
00200         }
00201     }
00202     else
00203     {
00204         // The general case
00205         // a != 0; b != 0; c != 0
00206         testpoint = -(b*ybmin + c*zbmin + d)/a;
00207         if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
00208         {
00209             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
00210         }
00211         if ( intersection_points.size() == 6 ) return;
00212         testpoint = -(b*ybmax + c*zbmin + d)/a;
00213         if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
00214         {
00215             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
00216         }
00217         if ( intersection_points.size() == 6 ) return;
00218         testpoint = -(b*ybmax + c*zbmax + d)/a;
00219         if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
00220         {
00221             add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
00222         }
00223         if ( intersection_points.size() == 6 ) return;
00224         testpoint = -(b*ybmin + c*zbmax + d)/a;
00225         if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
00226         {
00227             add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
00228         }
00229         if ( intersection_points.size() == 6 ) return;
00230         testpoint = -(a*xbmin + c*zbmin + d)/b;
00231         if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
00232         {
00233             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
00234         }
00235         if ( intersection_points.size() == 6 ) return;
00236         testpoint = -(a*xbmax + c*zbmin + d)/b;
00237         if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
00238         {
00239             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
00240         }
00241         if ( intersection_points.size() == 6 ) return;
00242         testpoint = -(a*xbmax + c*zbmax + d)/b;
00243         if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
00244         {
00245             add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
00246         }
00247         if ( intersection_points.size() == 6 ) return;
00248         testpoint = -(a*xbmin + c*zbmax + d)/b;
00249         if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
00250         {
00251             add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
00252         }
00253         if ( intersection_points.size() == 6 ) return;
00254         testpoint = -(a*xbmin + b*ybmin + d)/c;
00255         if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
00256         {
00257             add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
00258         }
00259         if ( intersection_points.size() == 6 ) return;
00260         testpoint = -(a*xbmax + b*ybmin + d)/c;
00261         if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
00262         {
00263             add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
00264         }
00265         if ( intersection_points.size() == 6 ) return;
00266         testpoint = -(a*xbmax + b*ybmax + d)/c;
00267         if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
00268         {
00269             add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
00270         }
00271         if ( intersection_points.size() == 6 ) return;
00272         testpoint = -(a*xbmin + b*ybmax + d)/c;
00273         if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
00274         {
00275             add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
00276         }
00277     }
00278 }
00279 
00280 //===========================================================================
00281 //Function Name: add_unique_point
00282 //Description:   Given a point and a list of points, add the point to the list
00283 //               unless it is already in the list within GEOMETRY_RESABS.
00284 //Author:
00285 //Date: 3/2005
00286 //===========================================================================
00287 void FBDataUtil::add_unique_point
00288 (
00289     DLIList<CubitVector > &points,
00290     const CubitVector &pt
00291 )
00292 {
00293     int ipt;
00294     for ( ipt = 0; ipt < points.size(); ipt++ )
00295     {
00296         double dist = pt.distance_between( points[ipt] );
00297         if ( dist <= GEOMETRY_RESABS )
00298         {
00299             return;
00300         }
00301     }
00302     points.append( pt );
00303 }
00304 
00305 //===========================================================================
00306 //Function Name: FBmake_xy_plane
00307 //Description:   Makes a triangle plane centered on the origin in the x-y 
00308 //               direction.  
00309 //Author: jdfowle
00310 //Date: 1/2004
00311 //===========================================================================
00312 CubitStatus FBDataUtil::FBmake_xy_plane(std::vector<double>& verts, 
00313                                         std::vector<int>& conns, 
00314                                         double xsize, 
00315                                         double ysize, 
00316                                         int numx, 
00317                                         int numy)
00318 {
00319 int i, j;
00320 double xc, yc, zc, xinc, yinc;
00321 
00322   zc = 0.;
00323   xinc = xsize/(double)(numx-1);
00324   yinc = ysize/(double)(numy-1);
00325   xc = -0.5*xsize;
00326 //  Make the coordinates
00327   for ( i = 0; i < numx; i++ ) {
00328     yc = -0.5*ysize;
00329     for ( j = 0; j < numy; j++ ) {    
00330        verts.push_back(xc);
00331        verts.push_back(yc);
00332        verts.push_back(zc);
00333        yc += yinc;
00334     }
00335     xc += xinc; 
00336   }
00337 //  Make the connections
00338 int i1, i2, i3, i4;
00339 
00340   for ( i = 0; i < numx-1; i++ ) {
00341   
00342     for ( j = 0; j < numy-1; j++ ) {
00343       i1 = numy*i + j;
00344       i2 = numy*(i+1) + j;
00345       i3 = numy*i + j + 1;
00346       i4 = numy*(i+1) + j + 1;
00347       conns.push_back(i1);
00348       conns.push_back(i2);
00349       conns.push_back(i4);
00350       conns.push_back(i1);
00351       conns.push_back(i4);
00352       conns.push_back(i3);          
00353     }
00354   }
00355   
00356   return CUBIT_SUCCESS;
00357 }
00358 
00359 //===========================================================================
00360 //Function Name: rotate_FB_object
00361 //Description:   Rotates a facetedboolean object so that the z-direction 
00362 //               points in the direction of NormalDir
00363 //Author: jdfowle
00364 //Date: 1/2004
00365 //===========================================================================
00366 CubitStatus FBDataUtil::rotate_FB_object(std::vector<double>& verts,
00367                              CubitVector &NormalDir, CubitVector &CenterPt)
00368 {
00369 unsigned int i; 
00370 double l1, l2, l3, m1, m2, m3, n1, n2, n3;
00371 double tnx, tny, tnz, temp;
00372 double tx, ty, tz;
00373 double xpcen, ypcen, zpcen;
00374 
00375   tnx = NormalDir.x();
00376   tny = NormalDir.y();
00377   tnz = NormalDir.z();
00378   xpcen = CenterPt.x();
00379   ypcen = CenterPt.y();
00380   zpcen = CenterPt.z();
00381   
00382   n1 = tnx; n2 = tny; n3 = tnz;
00383   l1 = n3; l2 = 0.; l3 = -n1;
00384   m1 = -n1*n2; m2 = n1*n1 + n3*n3; m3 = -n2*n3;
00385   temp = sqrt(l1*l1 + l3*l3);
00386   if ( fabs(temp) > 1.e-6 ) {
00387     l1 /= temp; l3 /= temp;
00388     temp = sqrt(m1*m1 + m2*m2 + m3*m3);
00389     m1 /= temp; m2 /= temp; m3 /= temp;
00390   } else {
00391     l1 = 1.; l2 = l3 = 0.;
00392     m1 = m2 = 0.; m3 = 1.;
00393     n1 = n3 = 0.; n2 = -1.;
00394   }
00395 //  Rotate the plane, whiuch is normal to the Z axis and through the origin,
00396 //  so that it will be normal to threepointnormal.  Also translate it so
00397 // that the center point, (0,0,0), moves to threepointcenter.
00398   for ( i = 0; i < verts.size(); i += 3 ) {
00399     tx = verts[i];
00400     ty = verts[i+1];
00401     tz = verts[i+2];
00402     verts[i] = tx*l1 + ty*m1 + tz*n1 + xpcen;
00403     verts[i+1] = tx*l2 + ty*m2 + tz*n2 + ypcen;
00404     verts[i+2] = tx*l3 + ty*m3 + tz*n3 + zpcen;
00405   } 
00406   
00407   return CUBIT_SUCCESS;
00408 }
00409 
00410 int FBDataUtil::makeahashvaluefrom_coord(double x, double y, double z, int numhashbins)
00411 {
00412 double sum;
00413 
00414       if ( fabs(x) < 1.e-3 ) x = 0.0;
00415       if ( fabs(y) < 1.e-3 ) y = 0.0;
00416       if ( fabs(z) < 1.e-3 ) z = 0.0;
00417       sum = (int)(10000.0*fabs(x) + 0.5) + 
00418                     (int)(10000.0*fabs(y) + 0.5) + 
00419                     (int)(10000.0*fabs(z) + 0.5);
00420       
00421       return (int)(sum) % numhashbins;
00422 }
00423  
00424 CubitStatus FBDataUtil::FBmake_cylinder(std::vector<double>& verts, 
00425                                         std::vector<int>& coords, 
00426                                         double radius, 
00427                                         double length, 
00428                                         int nr, 
00429                                         int nz)
00430 {
00431   int i, j;
00432   CubitStatus status;   
00433   double cfac, rinc, linc;
00434   double x, y, z;
00435   int istart, iend, V3, pend;
00436   double zoffset, lpos, rpos, xrad, yrad;
00437   
00438   status = CUBIT_SUCCESS;
00439   
00440   rinc = 360.0/(double)nr;
00441   linc = length/(double)nz;
00442   cfac = CUBIT_PI/180.;
00443   
00444     istart = 0; iend = nz+1;
00445     V3 = (nz+1)*nr;
00446     pend = nz;
00447   
00448   
00449   //  Make the points.
00450   
00451   zoffset = 0.0;
00452   lpos = -0.5*length; 
00453   xrad = radius;
00454   yrad = radius;
00455   for ( i = istart; i < iend; i++ ) {
00456     rpos = 10.0;
00457     for ( j = 0; j < nr; j++ ) {
00458       x = xrad*cos(cfac*rpos);
00459       y = yrad*sin(cfac*rpos);
00460       z = lpos;
00461       verts.push_back(x);
00462       verts.push_back(y);
00463       verts.push_back(z);      
00464       rpos += rinc;   
00465     }
00466     lpos += linc;
00467     zoffset += linc;
00468   } 
00469   //  Add the two points on the axis at the ends.
00470   verts.push_back(0.);
00471   verts.push_back(0.);
00472   verts.push_back(-0.5*length);
00473   verts.push_back(0.);
00474   verts.push_back(0.);
00475   verts.push_back(0.5*length);
00476     
00477   //  Make the triangles.
00478   int vertnum;
00479   vertnum = 0;
00480   for ( i = 0; i < pend; i++ ) {
00481     for ( j = 0; j < nr-1; j++ ) {
00482 //      facet_ptr = new CubitFacetData( points[vertnum+j],points[vertnum+j+1], points[vertnum+j+nr] );
00483     coords.push_back(vertnum+j);
00484     coords.push_back(vertnum+j+1);
00485     coords.push_back(vertnum+j+nr);
00486     
00487 //      facet_ptr = new CubitFacetData( points[vertnum+j+1],points[vertnum+j+1+nr], points[vertnum+j+nr] );
00488     coords.push_back(vertnum+j+1);
00489     coords.push_back(vertnum+j+1+nr);
00490     coords.push_back(vertnum+j+nr);
00491  
00492     }
00493 //    facet_ptr = new CubitFacetData( points[vertnum],points[vertnum+nr], points[vertnum+2*nr-1] );
00494     coords.push_back(vertnum);
00495     coords.push_back(vertnum+nr);
00496     coords.push_back(vertnum+2*nr-1);
00497 
00498 //    facet_ptr = new CubitFacetData( points[vertnum+nr-1],points[vertnum], points[vertnum+2*nr-1] );
00499     coords.push_back(vertnum+nr-1);
00500     coords.push_back(vertnum);
00501     coords.push_back(vertnum+2*nr-1);
00502     
00503     vertnum += nr;
00504   }
00505   
00506   //  Endcap(s)
00507   for ( i = 0; i < nr-1; i++ ) { // top cap
00508 //    facet_ptr = new CubitFacetData( points[vertnum+i],points[vertnum+i+1], points[V3+1] );
00509     coords.push_back(vertnum+i);
00510     coords.push_back(vertnum+i+1);
00511     coords.push_back(V3+1);
00512 
00513   }   
00514 //  facet_ptr = new CubitFacetData( points[nr-1+vertnum],points[vertnum], points[V3+1] );
00515     coords.push_back(vertnum+nr-1);
00516     coords.push_back(vertnum);
00517     coords.push_back(V3+1);
00518   
00519   
00520   for ( i = 0; i < nr-1; i++ ) { // bottom cap
00521 //    facet_ptr = new CubitFacetData( points[i+1],points[i], points[V3] );
00522     coords.push_back(i+1);
00523     coords.push_back(i);
00524     coords.push_back(V3);
00525  
00526   }   
00527 //  facet_ptr = new CubitFacetData( points[0],points[nr-1], points[V3] );
00528     coords.push_back(0);
00529     coords.push_back(nr-1);
00530     coords.push_back(V3);
00531  
00532   return status;
00533   
00534 }
00535 
00536 double FBDataUtil::project_point_to_plane(double *point, double a, double b, 
00537                                           double c, double d, 
00538                                           double *projected_pt)
00539 {
00540 double signed_distance;
00541 
00542   signed_distance = point[0]*a + point[1]*b + point[2]*c + d;
00543   projected_pt[0] = point[0] - signed_distance*a;
00544   projected_pt[1] = point[1] - signed_distance*b;
00545   projected_pt[2] = point[2] - signed_distance*c;
00546 
00547 
00548   return signed_distance;
00549 
00550 }
00551 
00552 
00553 #define EPS_SEG_SEG 1.e-6
00554 //  From Schneider and Eberly,"Geometric Tools for Computer Graphics",
00555 //  Chap. 10.8, segment to segment
00556 //  Note:  if the segments are parallel, *s and *t are undefined; *parallel
00557 //  is true; *sunclipped  = 0.0, and *tunclipped is the corresponding value
00558 //  for t.
00559 double FBDataUtil::closest_seg_seg_dist(double *p0, double *d0, double *p1, 
00560                              double *d1, double *s, double *t, 
00561                              double *sunclipped, double *tunclipped, 
00562                              bool *parallel)
00563 {
00564 double ux, uy, uz;
00565 double a, b, c, d, e, det;
00566 double snum, sdenom, tnum, tdenom;
00567 
00568   *parallel = false;
00569   
00570   ux = p0[0] - p1[0];
00571   uy = p0[1] - p1[1];
00572   uz = p0[2] - p1[2];
00573   
00574   a = d0[0]*d0[0] + d0[1]*d0[1] + d0[2]*d0[2];
00575   b = d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2];
00576   c = d1[0]*d1[0] + d1[1]*d1[1] + d1[2]*d1[2];
00577   d = d0[0]*ux + d0[1]*uy + d0[2]*uz;  
00578   e = d1[0]*ux + d1[1]*uy + d1[2]*uz;  
00579   det = a*c - b*b;
00580 
00581   if ( det < EPS_SEG_SEG ) {
00582     snum = 0.;
00583     tnum = e;
00584     tdenom = c;
00585     sdenom = 1.;
00586     *sunclipped = snum/sdenom;
00587     *tunclipped = tnum/tdenom;    
00588     double f = ux*ux + uy*uy + uz*uz;
00589     *parallel = true; 
00590     return sqrt(f - e*e/c); 
00591   } else {
00592     snum = b*e - c*d;
00593     tnum = a*e - b*d;
00594     sdenom = det;
00595     *sunclipped = snum/det;
00596     *tunclipped = tnum/det;    
00597   }
00598 
00599 
00600   if ( snum < 0.0 ) {
00601     snum = 0.0;
00602     tnum = e;
00603     tdenom = c;
00604   } else if ( snum > det ) {
00605     snum = det;
00606     tnum = e + b;
00607     tdenom = c;
00608   } else {
00609     tdenom = det;
00610   }
00611   
00612   if ( tnum < 0.0 ) {
00613     tnum = 0.0;
00614     if ( -d < 0.0 ) {
00615       snum = 0.0;
00616     } else if ( -d > a ) {
00617         snum = sdenom;
00618     } else {
00619         snum = -d;
00620         sdenom = a;
00621     }
00622   } else if ( tnum > tdenom ) {
00623       tnum = tdenom;
00624     if ( (-d + b) < 0.0 ) {
00625         snum = 0.0;
00626     } else if ( (-d + b) > a ) {
00627           snum = sdenom;
00628     } else {
00629           snum = -d + b;
00630           sdenom = a;
00631     }
00632   }
00633 
00634   *s = snum/sdenom;
00635   *t = tnum/tdenom;
00636 
00637   double vx, vy, vz;
00638   vx = p0[0] + (*s*d0[0]) - (p1[0] + (*t*d1[0]));
00639   vy = p0[1] + (*s*d0[1]) - (p1[1] + (*t*d1[1]));
00640   vz = p0[2] + (*s*d0[2]) - (p1[2] + (*t*d1[2]));
00641 
00642   return sqrt(vx*vx + vy*vy + vz*vz);
00643        
00644 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines