cgma
Cholla.cpp
Go to the documentation of this file.
00001 //- Class:       Cholla
00002 //- Description: C Interface for the Cholla module
00003 //-              Facet-based geometry definition and evaluation
00004 //- Owner:       Steven J. Owen
00005 //- Checked by:
00006 //- Version:
00007 
00008 #include "time.h"
00009 #include "Cholla.h"
00010 #include "ChollaEngine.hpp"
00011 #include "CubitPoint.hpp"
00012 #include "CubitPointData.hpp"
00013 #include "CubitFacetEdge.hpp"
00014 #include "CubitFacetEdgeData.hpp"
00015 #include "CubitFacet.hpp"
00016 #include "CubitFacetData.hpp"
00017 #include "CubitQuadFacet.hpp"
00018 #include "CubitQuadFacetData.hpp"
00019 #include "CastTo.hpp"
00020 #include "TDFacetBoundaryPoint.hpp"
00021 #include "FacetEvalTool.hpp"
00022 #include "FacetDataUtil.hpp"
00023 
00024 #include <iostream>
00025 #include <fstream>
00026 using std::ifstream;
00027 
00028 #define END_OF_BLOCKS "END_OF_BLOCKS"
00029 
00030 // Internal Static function prototypes
00031 
00032 static void time_stamp( FILE *fp );
00033 static CubitStatus build_facets(int numTri, int numQuad, int numEdge, 
00034                                 int numVert, int* triEdge, 
00035                                 int *quadEdge, 
00036                                 int* edgeVert, double* vert,
00037                                 CubitPoint **point_array, 
00038                                 CubitFacetEdge **edge_array, 
00039                                 CubitFacet **facet_array,
00040                                 CubitQuadFacet **qfacet_array,
00041                                 DLIList <FacetEntity *> &facet_list);
00042 static CubitStatus extractEdgeTangsFromFile( int num_file_edges,
00043                                              int *edge_nodes_file,
00044                                              double *edge_node_tangs_file,
00045                                              int numTri,
00046                                              int *edgeVert,
00047                                              double *edgeVertTang );
00048 static CubitStatus extractTriNormalsFromFile( int num_file_tris,
00049                                               int *tri_nodes_file,
00050                                               double *tri_node_norms_file,
00051                                               int numTri,
00052                                               int *edgeVert,
00053                                               int *triEdge,
00054                                               double *triVertNorms );
00055 static CubitStatus extractQuadNormalsFromFile( int num_file_quads,
00056                                                int *quad_nodes_file,
00057                                                double *quad_node_norms_file,
00058                                                int numQuad,
00059                                                int *edgeVert,
00060                                                int *quadEdge,
00061                                                double *quadVertNorms );
00062 static void constructTriVerts( int triEdge[3],
00063                                int *edgeVert,
00064                                int this_tri[3] );
00065 static void constructQuadVerts( int quadEdge[3],
00066                                 int *edgeVert,
00067                                 int this_quad[4] );
00068 static void checkMemoryAllocations( int num_nodes_per_elem,
00069                                     int additional_num_elems,
00070                                     int *num_elems,
00071                                     int **nodes,
00072                                     double **node_normals );
00073 
00074 static int classify_local_convexity_at_edge(CubitFacetEdge *edge_ptr);
00075 
00076 //===========================================================================
00077 //  Function: constructBezier
00078 //  Purpose:  Return control points for a quartic Bezier for each face, using
00079 //            a given feature angle tolerance to distinguish C0 and C1 edges.
00080 //            Separate arrays of tris and quad faces should be supplied
00081 //  Date:     10/28/2002
00082 //  Author:   sjowen
00083 //===========================================================================
00084 void constructBezier( double angle, int numTri, int numQuad, int numEdge,
00085                       int numVert, int* triEdge, int* quadEdge,
00086                       int* edgeVert, double* vert,
00087                       double* edgeCtrlPts, double* triCtrlPts,
00088                       double* quadCtrlPts )
00089 {
00090 
00091   // create arrays of facet entities from the arrays
00092 
00093   CubitPoint **point_array = new CubitPoint * [numVert];
00094   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
00095   CubitFacet **facet_array = NULL;
00096   if (numTri > 0)
00097     facet_array = new CubitFacet * [numTri];
00098   CubitQuadFacet **qfacet_array = NULL;
00099   if (numQuad > 0)
00100     qfacet_array = new CubitQuadFacet * [numQuad];
00101   DLIList<FacetEntity *> facet_list;
00102   CubitStatus stat = CUBIT_SUCCESS;
00103   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
00104                       edgeVert, vert, point_array, edge_array, facet_array,
00105                       qfacet_array, facet_list);
00106   if (CUBIT_SUCCESS != stat) {
00107     PRINT_ERROR("Failed to build facets.\n");
00108     return;
00109   }
00110 
00111   // create lists of points and edges
00112 
00113   int ii;
00114   DLIList<CubitPoint *> point_list;
00115   for (ii=0; ii<numVert; ii++)  
00116     point_list.append(point_array[ii]);  
00117 
00118   DLIList<CubitFacetEdge *> edge_list;
00119   for (ii=0; ii<numEdge; ii++)
00120     edge_list.append(edge_array[ii]);
00121 
00122   // prepare the ChollaEngine
00123 
00124   DLIList<FacetEntity *> fpoint_list, fedge_list;
00125   CAST_LIST_TO_PARENT( point_list, fpoint_list );
00126   CAST_LIST_TO_PARENT( edge_list,  fedge_list );
00127   ChollaEngine *ceng = new ChollaEngine( facet_list, fedge_list, fpoint_list );
00128 
00129   CubitBoolean use_feature_angle;
00130   if (angle < 0.00001 || angle > 179.9999)
00131     use_feature_angle = CUBIT_FALSE;
00132   else
00133     use_feature_angle = CUBIT_TRUE;
00134   int interp_order = 4;  // assume we are creating beziers
00135   CubitBoolean smooth_non_manifold = CUBIT_TRUE;
00136   CubitBoolean split_surfaces = CUBIT_FALSE;
00137 
00138   // create the geometry
00139     
00140   //CubitStatus rv = 
00141   ceng->create_geometry(use_feature_angle,angle,
00142                         interp_order,smooth_non_manifold,
00143                         split_surfaces);
00144 
00145   // get the control points on edges
00146 
00147   int jj;
00148   int index;
00149   CubitVector *control_points;
00150   CubitFacetEdge *edge_ptr;
00151   for (ii=0; ii<numEdge; ii++)
00152   {
00153     edge_ptr = edge_array[ii];
00154     control_points = edge_ptr->control_points();
00155     assert(control_points != NULL);
00156     for (jj=0; jj<NUM_EDGE_CPTS; jj++)
00157     {
00158       index = 3 * (ii * NUM_EDGE_CPTS + jj); 
00159       edgeCtrlPts[index] = control_points[jj].x();
00160       edgeCtrlPts[index+1] = control_points[jj].y();
00161       edgeCtrlPts[index+2] = control_points[jj].z();
00162     }
00163   }
00164 
00165   // get the control points on triangles
00166 
00167   CubitFacet *facet_ptr;
00168   for (ii=0; ii<numTri; ii++)
00169   {
00170     facet_ptr = facet_array[ii];
00171     control_points = facet_ptr->control_points();
00172     assert(control_points != NULL);
00173     for (jj=0; jj<NUM_TRI_CPTS; jj++)
00174     {
00175       index = 3 * (ii * NUM_TRI_CPTS + jj); 
00176       triCtrlPts[index] = control_points[jj].x();
00177       triCtrlPts[index+1] = control_points[jj].y();
00178       triCtrlPts[index+2] = control_points[jj].z();
00179     }
00180   }
00181 
00182     // get the control points on quads
00183 
00184   //CubitVector qctrl_points[NUM_QUAD_CPTS];
00185   CubitQuadFacet *qfacet_ptr;
00186   for (ii=0; ii<numQuad; ii++)
00187   {
00188     qfacet_ptr = qfacet_array[ii];
00189     index = 3 * (ii * NUM_QUAD_CPTS);
00190     qfacet_ptr->get_control_points( &(quadCtrlPts[index]) );
00191   }
00192 
00193   // delete the ChollaEngine
00194 
00195   ceng->delete_eval_tools();
00196   ceng->delete_me();
00197   delete ceng;
00198 
00199   // delete temp arrays (facets are deleted with the eval tool)
00200 
00201   for (ii=0; ii<numQuad; ii++)
00202   {
00203     qfacet_array[ii]->remove_tri_facets();
00204     delete qfacet_array[ii];
00205   }
00206   delete [] point_array;
00207   delete [] edge_array;
00208   if (facet_array)
00209     delete [] facet_array; 
00210   if (qfacet_array)
00211     delete [] qfacet_array; 
00212 
00213 }
00214 
00215 //===========================================================================
00216 //  Function: constructTriNormals
00217 //  Purpose:  Return normals and tangents for a set of triangles using
00218 //            a given feature angle tolerance to distinguish C0 and C1 edges.
00219 //  Date:     10/10/2003
00220 //  Author:   sjowen
00221 //===========================================================================
00222 void constructTriNormals( double angle, int numTri, int numEdge,
00223                           int numVert, int* triEdge,
00224                           int* edgeVert, double* vert,
00225                           double* triVertNorms, double* edgeVertTang,
00226                           int* edgeTypeFlags)
00227 {
00228 
00229   // create arrays of facet entities from the arrays
00230 
00231   int* quadEdge = NULL;
00232   CubitPoint **point_array = new CubitPoint * [numVert];
00233   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
00234   CubitFacet **facet_array = NULL;
00235   if (numTri > 0)
00236     facet_array = new CubitFacet * [numTri];
00237   int numQuad = 0;
00238   CubitQuadFacet **qfacet_array = NULL;
00239   DLIList<FacetEntity *> facet_list;
00240   CubitStatus stat = CUBIT_SUCCESS;
00241   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
00242                       edgeVert, vert, point_array, edge_array, facet_array,
00243                       qfacet_array, facet_list);
00244   if (CUBIT_SUCCESS != stat) {
00245     PRINT_ERROR("Failed to build facets.\n");
00246     return;
00247   }
00248 
00249   // create lists of points and edges
00250 
00251   int ii;
00252   DLIList<CubitPoint *> point_list;
00253   for (ii=0; ii<numVert; ii++)  
00254     point_list.append(point_array[ii]);  
00255 
00256   DLIList<CubitFacetEdge *> edge_list;
00257   for (ii=0; ii<numEdge; ii++){
00258     edge_list.append(edge_array[ii]);
00259   }
00260   
00261 
00262   // prepare the ChollaEngine
00263 
00264   DLIList<FacetEntity *> fpoint_list, fedge_list;
00265   CAST_LIST_TO_PARENT( point_list, fpoint_list );
00266   CAST_LIST_TO_PARENT( edge_list,  fedge_list );
00267   ChollaEngine *ceng = new ChollaEngine( facet_list, fedge_list, fpoint_list );
00268 
00269   CubitBoolean use_feature_angle;
00270   if (angle < 0.00001 || angle > 179.9999)
00271     use_feature_angle = CUBIT_FALSE;
00272   else
00273     use_feature_angle = CUBIT_TRUE;
00274   int interp_order = 4;  // assume we are creating beziers
00275   CubitBoolean smooth_non_manifold = CUBIT_TRUE;
00276   CubitBoolean split_surfaces = CUBIT_FALSE;
00277 
00278   // create the geometry
00279     
00280   //CubitStatus rv = 
00281   ceng->create_geometry(use_feature_angle,angle,
00282                         interp_order,smooth_non_manifold,
00283                         split_surfaces);
00284 
00285   // get the control points on edges and infer the tangents
00286 
00287   int index = 0;
00288   CubitVector *control_points;
00289   CubitFacetEdge *edge_ptr;
00290   CubitVector v0, v1, t0, t1;
00291   for (ii=0; ii<numEdge; ii++)
00292   {
00293     edge_ptr = edge_array[ii];
00294     v0 = edge_ptr->point(0)->coordinates();
00295     v1 = edge_ptr->point(1)->coordinates();
00296     control_points = edge_ptr->control_points();
00297     assert(control_points != NULL);
00298     index = 6 * ii;
00299     t0 = control_points[0] - v0;
00300     t1 = v1 - control_points[2];
00301     if(edge_ptr->is_flipped()){
00302       t0 *= (-1.0);
00303       t1 *= (-1.0);
00304     }
00305     t0.normalize();
00306     t1.normalize();
00307     edgeVertTang[index  ] = t0.x();
00308     edgeVertTang[index+1] = t0.y();
00309     edgeVertTang[index+2] = t0.z();
00310     edgeVertTang[index+3] = t1.x();
00311     edgeVertTang[index+4] = t1.y();
00312     edgeVertTang[index+5] = t1.z();
00313     edgeTypeFlags[ii]=classify_local_convexity_at_edge(edge_ptr);
00314   }
00315 
00316   // get the normal on triangles
00317 
00318   CubitVector n0, n1, n2;
00319   CubitFacet *facet_ptr;
00320   for (ii=0; ii<numTri; ii++)
00321   {
00322     facet_ptr = facet_array[ii];
00323     index = 9 * ii;
00324     n0 = facet_ptr->point(0)->normal(facet_ptr);
00325     n1 = facet_ptr->point(1)->normal(facet_ptr);
00326     n2 = facet_ptr->point(2)->normal(facet_ptr);
00327     triVertNorms[index    ] = n0.x();
00328     triVertNorms[index + 1] = n0.y();
00329     triVertNorms[index + 2] = n0.z();
00330     triVertNorms[index + 3] = n1.x();
00331     triVertNorms[index + 4] = n1.y();
00332     triVertNorms[index + 5] = n1.z();
00333     triVertNorms[index + 6] = n2.x();
00334     triVertNorms[index + 7] = n2.y();
00335     triVertNorms[index + 8] = n2.z();
00336   }
00337 
00338   // delete the ChollaEngine
00339 
00340   ceng->delete_eval_tools();
00341   ceng->delete_me();
00342   delete ceng;
00343 
00344   // delete temp arrays (facets are deleted with the eval tool)
00345 
00346   for (ii=0; ii<numQuad; ii++)
00347   {
00348     qfacet_array[ii]->remove_tri_facets();
00349     delete qfacet_array[ii];
00350   }
00351   delete [] point_array;
00352   delete [] edge_array;
00353   if (facet_array)
00354     delete [] facet_array; 
00355   if (qfacet_array)
00356     delete [] qfacet_array; 
00357 }
00358 
00359 //===========================================================================
00360 //  Function: constructQuadNormals
00361 //  Purpose:  Return normals and tangents for a set of quads using
00362 //            a given feature angle tolerance to distinguish C0 and C1 edges.
00363 //  Date:     10/10/2003
00364 //  Author:   sjowen
00365 //===========================================================================
00366 void constructQuadNormals( double angle, int numQuad, int numEdge,
00367                            int numVert, int* quadEdge,
00368                            int* edgeVert, double* vert,
00369                            double* quadVertNorms, double* edgeVertTang,
00370                            int* edgeTypeFlags )
00371 {
00372 
00373   // create arrays of facet entities from the arrays
00374 
00375   int *triEdge = NULL;
00376   CubitPoint **point_array = new CubitPoint * [numVert];
00377   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
00378   CubitFacet **facet_array = NULL;
00379   CubitQuadFacet **qfacet_array = NULL; 
00380   int numTri = 0;
00381   if (numQuad > 0)
00382     qfacet_array = new CubitQuadFacet * [numQuad];
00383   DLIList<FacetEntity *> facet_list;
00384   CubitStatus stat = CUBIT_SUCCESS;
00385   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
00386                       edgeVert, vert, point_array, edge_array, facet_array,
00387                       qfacet_array, facet_list);
00388   if (CUBIT_SUCCESS != stat) {
00389     PRINT_ERROR("Failed to build facets.\n");
00390     return;
00391   }
00392   // create lists of points and edges
00393 
00394   int ii;
00395   
00396   DLIList<CubitPoint *> point_list;
00397   for (ii=0; ii<numVert; ii++)  
00398     point_list.append(point_array[ii]);  
00399 
00400   DLIList<CubitFacetEdge *> edge_list;
00401   for (ii=0; ii<numEdge; ii++)
00402     edge_list.append(edge_array[ii]);
00403 
00404   // prepare the ChollaEngine
00405 
00406   DLIList<FacetEntity *> fpoint_list, fedge_list;
00407   CAST_LIST_TO_PARENT( point_list, fpoint_list );
00408   CAST_LIST_TO_PARENT( edge_list,  fedge_list );
00409   ChollaEngine *ceng = new ChollaEngine( facet_list, fedge_list, fpoint_list );
00410 
00411   CubitBoolean use_feature_angle;
00412   if (angle < 0.00001 || angle > 179.9999)
00413     use_feature_angle = CUBIT_FALSE;
00414   else
00415     use_feature_angle = CUBIT_TRUE;
00416   int interp_order = 4;  // assume we are creating beziers
00417   CubitBoolean smooth_non_manifold = CUBIT_TRUE;
00418   CubitBoolean split_surfaces = CUBIT_FALSE;
00419 
00420   // create the geometry
00421  
00422   ceng->create_geometry(use_feature_angle,angle,
00423                         interp_order,smooth_non_manifold,
00424                         split_surfaces);
00425 
00426   // get the control points on edges and infer the tangents
00427 
00428   int mydebug = 0;
00429   FILE *fp = NULL;
00430   if (mydebug)
00431   {
00432     fp = fopen("edges.txt", "w");
00433   }
00434   int index = 0;
00435   CubitVector *control_points;
00436   CubitFacetEdge *edge_ptr;
00437   CubitVector v0, v1, t0, t1;
00438   for (ii=0; ii<numEdge; ii++)
00439   {
00440     edge_ptr = edge_array[ii];
00441     v0 = edge_ptr->point(0)->coordinates();
00442     v1 = edge_ptr->point(1)->coordinates();
00443     control_points = edge_ptr->control_points();
00444     assert(control_points != NULL);
00445     if (mydebug)
00446     {
00447       int kk;
00448       for(kk=0; kk<3; kk++);
00449       fprintf(fp, "(%d) %.6f %.6f %.6f\n", edge_ptr->id(), v0.x(), v0.y(), v0.z());
00450       fprintf(fp, "    %.6f %.6f %.6f\n", control_points[0].x(), control_points[0].y(), control_points[0].z());
00451       fprintf(fp, "    %.6f %.6f %.6f\n", control_points[1].x(), control_points[1].y(), control_points[1].z());
00452       fprintf(fp, "    %.6f %.6f %.6f\n", control_points[2].x(), control_points[2].y(), control_points[2].z());
00453       fprintf(fp, "    %.6f %.6f %.6f\n", v1.x(), v1.y(), v1.z());
00454     }
00455     index = 6 * ii; 
00456     t0 = control_points[0] - v0;
00457     t1 = v1 - control_points[2];
00458     if(edge_ptr->is_flipped()){
00459       t0 *= (-1.0);
00460       t1 *= (-1.0);
00461     }
00462     t0.normalize();
00463     t1.normalize();
00464     edgeVertTang[index  ] = t0.x();
00465     edgeVertTang[index+1] = t0.y();
00466     edgeVertTang[index+2] = t0.z();
00467     edgeVertTang[index+3] = t1.x();
00468     edgeVertTang[index+4] = t1.y();
00469     edgeVertTang[index+5] = t1.z();
00470     edgeTypeFlags[ii]=classify_local_convexity_at_edge(edge_ptr);
00471   }
00472   if (mydebug)
00473   {
00474     fclose(fp);
00475   }
00476 
00477   // get the normal on quad
00478 
00479   CubitVector n0, n1, n2, n3;
00480   CubitQuadFacet *qfacet_ptr;
00481   for (ii=0; ii<numQuad; ii++)
00482   {
00483     qfacet_ptr = qfacet_array[ii];
00484     index = 12 * ii;
00485     n0 = qfacet_ptr->point(0)->normal(qfacet_ptr);
00486     n1 = qfacet_ptr->point(1)->normal(qfacet_ptr);
00487     n2 = qfacet_ptr->point(2)->normal(qfacet_ptr);
00488     n3 = qfacet_ptr->point(3)->normal(qfacet_ptr);
00489     quadVertNorms[index     ] = n0.x();
00490     quadVertNorms[index + 1 ] = n0.y();
00491     quadVertNorms[index + 2 ] = n0.z();
00492     quadVertNorms[index + 3 ] = n1.x();
00493     quadVertNorms[index + 4 ] = n1.y();
00494     quadVertNorms[index + 5 ] = n1.z();
00495     quadVertNorms[index + 6 ] = n2.x();
00496     quadVertNorms[index + 7 ] = n2.y();
00497     quadVertNorms[index + 8 ] = n2.z();
00498     quadVertNorms[index + 9 ] = n3.x();
00499     quadVertNorms[index + 10] = n3.y();
00500     quadVertNorms[index + 11] = n3.z();
00501   }
00502   
00503   // delete the ChollaEngine
00504 
00505   ceng->delete_eval_tools();
00506   ceng->delete_me();
00507   delete ceng;
00508 
00509   // delete temp arrays
00510 
00511   for (ii=0; ii<numQuad; ii++)
00512   {
00513     qfacet_array[ii]->remove_tri_facets();
00514     delete qfacet_array[ii];
00515   }
00516   delete [] point_array;
00517   delete [] edge_array;
00518   if (facet_array)
00519     delete [] facet_array; 
00520   if (qfacet_array)
00521     delete [] qfacet_array; 
00522 }
00523 
00524 //===========================================================================
00525 //  Function: build_facets
00526 //  Purpose:  static function that creates facet entities
00527 //            from initial the arrays
00528 //  Date:     10/28/2002
00529 //  Author:   sjowen
00530 //===========================================================================
00531 static CubitStatus build_facets(int numTri, int numQuad, int numEdge, 
00532                                 int numVert, int* triEdge, int *quadEdge, 
00533                                 int* edgeVert, double* vert,
00534                                 CubitPoint **point_array, 
00535                                 CubitFacetEdge **edge_array, 
00536                                 CubitFacet **facet_array,
00537                                 CubitQuadFacet **qfacet_array,
00538                                 DLIList <FacetEntity *> &facet_list)
00539 {
00540   // create the points
00541 
00542   double x, y, z;
00543   int ii;
00544   
00545   for (ii=0; ii<numVert; ii++)
00546   {
00547     x = vert[ii*3];
00548     y = vert[ii*3 + 1];
00549     z = vert[ii*3 + 2];
00550     point_array[ii] = (CubitPoint *) new CubitPointData( x, y, z );
00551     point_array[ii]->set_id(ii);
00552   }
00553 
00554   // create the edges
00555 
00556   int ip, jp;
00557   for (ii=0; ii<numEdge; ii++)
00558   {
00559     ip = edgeVert[ii*2];
00560     jp = edgeVert[ii*2 + 1];
00561     assert(ip < numVert && jp < numVert && ip >= 0 && jp >= 0);
00562     edge_array[ii] = (CubitFacetEdge *) new CubitFacetEdgeData( point_array[ip],
00563                                                                 point_array[jp] );
00564     edge_array[ii]->set_id(ii);
00565   }
00566 
00567   // create the tri facets
00568 
00569   int begin_face;
00570   int jj, iedge;
00571   CubitFacetEdge *edges[4];
00572   CubitFacet *facet_ptr = NULL;
00573 
00574   for (ii=0; ii<numTri; ii++)
00575   {
00576     begin_face = 3 * ii;
00577     for(jj=0; jj<3; jj++)
00578     {
00579       iedge = triEdge[begin_face+jj];
00580       edges[jj] = edge_array[iedge];
00581     }
00582     facet_ptr = (CubitFacet *)
00583       new CubitFacetData(edges[0], edges[1], edges[2]);
00584     facet_list.append(facet_ptr);
00585     facet_array[ii] = facet_ptr;
00586   }
00587 
00588   // create the quad facets
00589 
00590   CubitQuadFacet *qfacet_ptr = NULL;
00591   for (ii=0; ii<numQuad; ii++)
00592   {
00593     begin_face = 4 * ii;
00594     for(jj=0; jj<4; jj++)
00595     {
00596       iedge = quadEdge[begin_face+jj];
00597       edges[jj] = edge_array[iedge];
00598     }
00599     qfacet_ptr = (CubitQuadFacet *)
00600       new CubitQuadFacetData(edges[0], edges[1], edges[2], edges[3]);
00601     facet_list.append(qfacet_ptr->get_tri_facet(0));
00602     facet_list.append(qfacet_ptr->get_tri_facet(1));
00603     qfacet_array[ii] = qfacet_ptr;
00604   }
00605 
00606   int mydebug = 0;
00607   if (mydebug)
00608   {
00609     DLIList<CubitFacet*> myfacet_list;
00610     CAST_LIST(facet_list, myfacet_list, CubitFacet);
00611     FacetDataUtil::write_facets("C:\\CUBIT\\cholla_apps\\examples\\sierra\\test.facets",myfacet_list);
00612   }
00613 
00614   return CUBIT_SUCCESS;
00615 
00616 }
00617 
00618 //===========================================================================
00619 //  Function: evalBezierFace
00620 //  Purpose:  Evaluate a set of quartic Bezier patches at a specified 
00621 //            parametric location given its control points. Evaluates 
00622 //            location, normal and/or derivative. Expects list of either 
00623 //            quads or tris (not both)
00624 //
00625 //  Date:     10/28/2002
00626 //  Author:   sjowen
00627 //===========================================================================
00628 void evalBezierFace( int numFace, int numEdge, int numVert, int numVertPerFace, 
00629                      int* faceEdge, int* edgeVert, 
00630                      double* vert,  double* faceCtrlPts,
00631                      double* edgeCtrlPts, 
00632                      int numLocs, double* paramLocation, 
00633                      double* location, double* normal, double* deriv )
00634 {
00635     // create arrays of facet entities from the arrays
00636 
00637   int numTri = 0;
00638   int numQuad = 0;
00639   CubitPoint **point_array = new CubitPoint * [numVert];
00640   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
00641   CubitFacet **facet_array = NULL;
00642   int *tri_edge = (numVertPerFace == 3) ? faceEdge : NULL;
00643   int *quad_edge = (numVertPerFace == 4) ? faceEdge : NULL;
00644   if (numVertPerFace == 3)
00645   {
00646     facet_array = new CubitFacet * [numFace];
00647     numTri = numFace;
00648   }
00649   CubitQuadFacet **qfacet_array = NULL;
00650   if (numVertPerFace == 4)
00651   {
00652     qfacet_array = new CubitQuadFacet * [numFace];
00653     numQuad = numFace;
00654   }
00655   DLIList<FacetEntity *> facet_list;
00656   CubitStatus stat = CUBIT_SUCCESS;
00657   stat = build_facets(numTri, numQuad, numEdge, numVert, tri_edge, quad_edge, 
00658                       edgeVert, vert, point_array, edge_array, facet_array,
00659                       qfacet_array, facet_list);
00660   if (CUBIT_SUCCESS != stat) {
00661     PRINT_ERROR("Failed to build facets.\n");
00662     return;
00663   }
00664 
00665   // set the control points on edges
00666 
00667   int index, ii, jj;
00668   CubitFacetEdge *edge_ptr;
00669   for (ii=0; ii<numEdge; ii++)
00670   {
00671     edge_ptr = edge_array[ii];
00672     index = 3 * (ii * NUM_EDGE_CPTS);
00673     edge_ptr->set_control_points( &(edgeCtrlPts[index]) );
00674   }
00675 
00676   // set the control points on triangles
00677 
00678   CubitFacet *facet_ptr;
00679   for (ii=0; ii<numTri; ii++)
00680   {
00681     facet_ptr = facet_array[ii];
00682     index = 3 * (ii * NUM_TRI_CPTS); 
00683     facet_ptr->set_control_points( &(faceCtrlPts[index]) );
00684   }
00685 
00686     // set the control points on quads
00687 
00688   CubitQuadFacet *qfacet_ptr;
00689   for (ii=0; ii<numQuad; ii++)
00690   {
00691     qfacet_ptr = qfacet_array[ii];
00692     index = 3 * (ii * NUM_QUAD_CPTS);
00693     qfacet_ptr->set_control_points( &(faceCtrlPts[index]) );
00694   }
00695 
00696   // Do the evaluations
00697 
00698   for (ii=0; ii<numTri; ii++)
00699   {
00700     for (jj=0; jj<numLocs; jj++)
00701     {
00702       facet_ptr = facet_array[ii];
00703       CubitVector areacoord(paramLocation[jj*3],
00704                             paramLocation[jj*3+1],
00705                             paramLocation[jj*3+2]);
00706       CubitVector eval_point;
00707       CubitVector *eval_point_ptr = &eval_point;
00708       CubitVector eval_normal;
00709       CubitVector *eval_normal_ptr = NULL;
00710       if (normal != NULL)
00711       {
00712         eval_normal_ptr = &eval_normal;
00713       }
00714       facet_ptr->evaluate( areacoord, eval_point_ptr, eval_normal_ptr );
00715       index = 3 * ii * jj;
00716       location[index] = eval_point.x();
00717       location[index+1] = eval_point.y();
00718       location[index+2] = eval_point.z();
00719       if (normal != NULL)
00720       {
00721         normal[index] = eval_normal.x();
00722         normal[index+1] = eval_normal.y();
00723         normal[index+2] = eval_normal.z();
00724       }
00725     }
00726   }
00727 
00728   for (ii=0; ii<numQuad; ii++)
00729   {
00730     for (jj=0; jj<numLocs; jj++)
00731     {
00732       qfacet_ptr = qfacet_array[ii];
00733       CubitVector eval_point;
00734       CubitVector *eval_point_ptr = &eval_point;
00735       CubitVector eval_normal;
00736       CubitVector *eval_normal_ptr = NULL;
00737       if (normal != NULL)
00738         eval_normal_ptr = &eval_normal;
00739       CubitVector eval_du, eval_dv;
00740       CubitVector *eval_du_ptr = NULL;
00741       CubitVector *eval_dv_ptr = NULL;
00742       if (deriv != NULL)
00743       {
00744         eval_du_ptr = &eval_du;
00745         eval_dv_ptr = &eval_dv;
00746       }
00747 
00748       qfacet_ptr->evaluate( paramLocation[2*jj], paramLocation[2*jj+1],
00749                             eval_point_ptr, eval_normal_ptr,
00750                             eval_du_ptr, eval_dv_ptr);
00751       index = 3 * ii * jj;
00752       location[index] = eval_point.x();
00753       location[index+1] = eval_point.y();
00754       location[index+2] = eval_point.z();
00755       if (normal != NULL)
00756       {
00757         normal[index] = eval_normal.x();
00758         normal[index+1] = eval_normal.y();
00759         normal[index+2] = eval_normal.z();
00760       }
00761       if (deriv != NULL)
00762       {
00763         index = 6 * ii * jj;
00764         deriv[index] = eval_du.x();
00765         deriv[index+1] = eval_du.y();
00766         deriv[index+2] = eval_du.z();
00767         deriv[index+3] = eval_dv.x();
00768         deriv[index+4] = eval_dv.y();
00769         deriv[index+5] = eval_dv.z();
00770       }
00771     }
00772   }
00773   
00774   // delete the temp arrays
00775 
00776   for (ii=0; ii<numQuad; ii++)
00777     delete qfacet_array[ii];
00778   for (ii=0; ii<numTri; ii++)
00779     delete facet_array[ii];
00780   for (ii=0; ii<numEdge; ii++)
00781     delete edge_array[ii];
00782   for (ii=0; ii<numVert; ii++)
00783     delete point_array[ii];
00784   delete [] point_array;
00785   delete [] edge_array;
00786   if (facet_array)
00787     delete [] facet_array; 
00788   if (qfacet_array)
00789     delete [] qfacet_array; 
00790 }
00791 
00792 //===========================================================================
00793 //  Function: evalBezierFaceFromNorms
00794 //  Purpose:  This is the same as the evalBezierFace except it differes by 
00795 //            the input arguments.  This function takes a list of normals 
00796 //            and tangents at the face vertices and computes bezier control 
00797 //            points internally.  Normals and tangents should have been 
00798 //            computed previously in constructTriNormals or constructQuadNormals.  
00799 //            This function is not as computationally efficient as evalBezierFace 
00800 //            since it requires Bezier control points to be computed
00801 //            as part of the call - however, it is more memory efficient, 
00802 //            requiring fewer variables to be stored with the calling application. 
00803 //            The following argument list describes only those arguments that 
00804 //            differ from evalBezierFace above.
00805 //
00806 //  Date:     10/15/2003
00807 //  Author:   sjowen
00808 //===========================================================================
00809 void evalBezierFaceFromNorms( int numFace, int numEdge, int numVert, 
00810                               int numVertPerFace, int* faceEdge, int* edgeVert, 
00811                               double* vert, double* vertNorms, double* edgeVertTang,
00812                               int numLocs, double* paramLocation, 
00813                               double* location, double* normal, double* deriv )
00814 {
00815     // create arrays of facet entities from the arrays
00816 
00817   int numTri = 0;
00818   int numQuad = 0;
00819   CubitPoint **point_array = new CubitPoint * [numVert];
00820   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
00821   CubitFacet **facet_array = NULL;
00822   int *triEdge = (numVertPerFace == 3) ? faceEdge : NULL;
00823   int *quadEdge = (numVertPerFace == 4) ? faceEdge : NULL;
00824   if (numVertPerFace == 3)
00825   {
00826     facet_array = new CubitFacet * [numFace];
00827     numTri = numFace;
00828   }
00829   CubitQuadFacet **qfacet_array = NULL;
00830   if (numVertPerFace == 4)
00831   {
00832     qfacet_array = new CubitQuadFacet * [numFace];
00833     numQuad = numFace;
00834   }
00835   DLIList<FacetEntity *> facet_list;
00836   CubitStatus stat = CUBIT_SUCCESS;
00837   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
00838                       edgeVert, vert, point_array, edge_array, facet_array,
00839                       qfacet_array, facet_list);
00840   if (CUBIT_SUCCESS != stat) {
00841     PRINT_ERROR("Failed to build facets.\n");
00842     return;
00843   }
00844 
00845   // set the normals on triangle vertices
00846 
00847   int index, ii, jj;
00848   CubitPoint *pt;
00849   CubitFacet *facet_ptr;
00850   CubitVector pt_normal;
00851   for (ii=0; ii<numTri; ii++)
00852   {
00853     facet_ptr = facet_array[ii];
00854     for(jj=0; jj<3; jj++)
00855     {
00856       index = ii * 9 + (3 * jj);
00857       pt_normal.x( vertNorms[ index ] );
00858       pt_normal.y( vertNorms[ index + 1 ] );
00859       pt_normal.z( vertNorms[ index + 2 ] );
00860       pt = facet_ptr->point(jj);
00861       TDFacetBoundaryPoint::add_facet_boundary_point( pt, facet_ptr, pt_normal );
00862     }
00863   }
00864 
00865   // set the normals on quad vertices
00866 
00867   CubitQuadFacet *qfacet_ptr;
00868   for (ii=0; ii<numQuad; ii++)
00869   {
00870     qfacet_ptr = qfacet_array[ii];
00871     for(jj=0; jj<4; jj++)
00872     {
00873       index = ii * 12 + (3 * jj);
00874       pt_normal.x( vertNorms[ index ] );
00875       pt_normal.y( vertNorms[ index + 1 ] );
00876       pt_normal.z( vertNorms[ index + 2 ] );
00877       pt = qfacet_ptr->point(jj);
00878       TDFacetBoundaryPoint::add_facet_boundary_point( pt, qfacet_ptr, pt_normal );
00879     }
00880   }
00881 
00882   // set the control points on edges
00883 
00884   int mydebug = 0;
00885   FILE *fp = NULL;
00886   if (mydebug)
00887     fp = fopen("edges.txt", "w");
00888   CubitFacetEdge *edge_ptr;
00889   CubitVector N0, N1, T0, T1, P0, P1;
00890   CubitVector Pi[3];
00891   CubitPoint *pt0, *pt1;
00892   for (ii=0; ii<numEdge; ii++)
00893   {
00894     index = 6 * ii;
00895     edge_ptr = edge_array[ii];
00896     pt0 = edge_ptr->point(0);
00897     pt1 = edge_ptr->point(1);
00898     P0 = pt0->coordinates();
00899     P1 = pt1->coordinates();
00900     N0 = pt0->normal( edge_ptr );
00901     N1 = pt1->normal( edge_ptr );
00902     T0.x( edgeVertTang[ index ] );
00903     T0.y( edgeVertTang[ index + 1] );
00904     T0.z( edgeVertTang[ index + 2] );
00905     T1.x( edgeVertTang[ index + 3] );
00906     T1.y( edgeVertTang[ index + 4] );
00907     T1.z( edgeVertTang[ index + 5] );
00908     FacetEvalTool::init_edge_control_points(P0, P1, N0, N1, T0, T1, Pi);
00909     if (mydebug)
00910     {
00911       int kk;
00912       for(kk=0; kk<3; kk++);
00913       fprintf(fp, "(%d) %.6f %.6f %.6f\n", edge_ptr->id(), P0.x(), P0.y(), P0.z());
00914       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[0].x(), Pi[0].y(), Pi[0].z());
00915       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[1].x(), Pi[1].y(), Pi[1].z());
00916       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[2].x(), Pi[2].y(), Pi[2].z());
00917       fprintf(fp, "    %.6f %.6f %.6f\n", P1.x(), P1.y(), P1.z());
00918     }
00919     edge_ptr->control_points( Pi, 4 );
00920   }
00921   if (mydebug)
00922     fclose(fp);
00923 
00924   // set the control points on triangles
00925 
00926   for (ii=0; ii<numTri; ii++)
00927   {
00928     facet_ptr = facet_array[ii];
00929     facet_ptr->init_patch();
00930   }
00931 
00932     // set the control points on quads
00933 
00934   for (ii=0; ii<numQuad; ii++)
00935   {
00936     qfacet_ptr = qfacet_array[ii];
00937     qfacet_ptr->init_patch();
00938   }
00939 
00940   // Do the evaluations
00941 
00942   index = 0;
00943   int nindex = 0;
00944   for (ii=0; ii<numTri; ii++)
00945   {
00946     facet_ptr = facet_array[ii];
00947     for (jj=0; jj<numLocs; jj++)
00948     {  
00949       CubitVector areacoord(paramLocation[jj*3],
00950                             paramLocation[jj*3+1],
00951                             paramLocation[jj*3+2]);
00952       CubitVector eval_point;
00953       CubitVector *eval_point_ptr = &eval_point;
00954       CubitVector eval_normal;
00955       CubitVector *eval_normal_ptr = NULL;
00956       if (normal != NULL)
00957       {
00958         eval_normal_ptr = &eval_normal;
00959       }
00960       facet_ptr->evaluate( areacoord, eval_point_ptr, eval_normal_ptr );
00961       location[index++] = eval_point.x();
00962       location[index++] = eval_point.y();
00963       location[index++] = eval_point.z();
00964       if (normal != NULL)
00965       {
00966         normal[nindex++] = eval_normal.x();
00967         normal[nindex++] = eval_normal.y();
00968         normal[nindex++] = eval_normal.z();
00969       }
00970     }
00971   }
00972 
00973   int dindex = 0;
00974   nindex = index = 0;
00975   for (ii=0; ii<numQuad; ii++)
00976   {
00977     for (jj=0; jj<numLocs; jj++)
00978     {
00979       qfacet_ptr = qfacet_array[ii];
00980       CubitVector eval_point;
00981       CubitVector *eval_point_ptr = &eval_point;
00982       CubitVector eval_normal;
00983       CubitVector *eval_normal_ptr = NULL;
00984       if (normal != NULL)
00985         eval_normal_ptr = &eval_normal;
00986       CubitVector eval_du, eval_dv;
00987       CubitVector *eval_du_ptr = NULL;
00988       CubitVector *eval_dv_ptr = NULL;
00989       if (deriv != NULL)
00990       {
00991         eval_du_ptr = &eval_du;
00992         eval_dv_ptr = &eval_dv;
00993       }
00994 
00995       qfacet_ptr->evaluate( paramLocation[2*jj], paramLocation[2*jj+1],
00996                             eval_point_ptr, eval_normal_ptr,
00997                             eval_du_ptr, eval_dv_ptr);
00998       location[index++] = eval_point.x();
00999       location[index++] = eval_point.y();
01000       location[index++] = eval_point.z();
01001       if (normal != NULL)
01002       {
01003         normal[nindex++] = eval_normal.x();
01004         normal[nindex++] = eval_normal.y();
01005         normal[nindex++] = eval_normal.z();
01006       }
01007       if (deriv != NULL)
01008       {
01009         index = 6 * ii * jj;
01010         deriv[dindex++] = eval_du.x();
01011         deriv[dindex++] = eval_du.y();
01012         deriv[dindex++] = eval_du.z();
01013         deriv[dindex++] = eval_dv.x();
01014         deriv[dindex++] = eval_dv.y();
01015         deriv[dindex++] = eval_dv.z();
01016       }
01017     }
01018   }
01019   
01020   // delete the temp arrays
01021 
01022   for (ii=0; ii<numQuad; ii++)
01023     delete qfacet_array[ii];
01024   for (ii=0; ii<numTri; ii++)
01025     delete facet_array[ii];
01026   for (ii=0; ii<numEdge; ii++)
01027     delete edge_array[ii];
01028   for (ii=0; ii<numVert; ii++)
01029     delete point_array[ii];
01030   delete [] point_array;
01031   delete [] edge_array;
01032   if (facet_array)
01033     delete [] facet_array; 
01034   if (qfacet_array)
01035     delete [] qfacet_array; 
01036 }
01037 
01038 //===========================================================================
01039 //  Function: projToBezierFace
01040 //  Purpose:  Project a set of x-y-z locations to Bezier patches.  Finds the 
01041 //            closest point on one of the patches.  Evaluates location, 
01042 //            normal and/or derivative.  Expects list of either quads or 
01043 //            tris (not both)
01044 //
01045 //  Date:     12/08/2003
01046 //  Author:   sjowen
01047 //===========================================================================
01048 void projToBezierFace( int numFace, int numEdge, int numVert, int numVertPerFace, 
01049                      int* faceEdge, int* edgeVert, 
01050                      double* vert,  double* faceCtrlPts,
01051                      double* edgeCtrlPts, 
01052                      int numLocs, double* xyz, 
01053                      int specify_tol, double converge_tol,
01054                      double* xyzOnFace, double* normal, double* deriv )
01055 {
01056     // create arrays of facet entities from the arrays
01057 
01058   int numTri = 0;
01059   int numQuad = 0;
01060   CubitPoint **point_array = new CubitPoint * [numVert];
01061   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
01062   CubitFacet **facet_array = NULL;
01063   int *tri_edge = (numVertPerFace == 3) ? faceEdge : NULL;
01064   int *quad_edge = (numVertPerFace == 4) ? faceEdge : NULL;
01065   if (numVertPerFace == 3)
01066   {
01067     facet_array = new CubitFacet * [numFace];
01068     numTri = numFace;
01069   }
01070   CubitQuadFacet **qfacet_array = NULL;
01071   if (numVertPerFace == 4)
01072   {
01073     qfacet_array = new CubitQuadFacet * [numFace];
01074     numQuad = numFace;
01075   }
01076   DLIList<FacetEntity *> facetentity_list;
01077   CubitStatus stat = CUBIT_SUCCESS;
01078   stat = build_facets(numTri, numQuad, numEdge, numVert, quad_edge, tri_edge, 
01079                       edgeVert, vert, point_array, edge_array, facet_array,
01080                       qfacet_array, facetentity_list);
01081   if (CUBIT_SUCCESS != stat) {
01082     PRINT_ERROR("Failed to build facets.\n");
01083     return;
01084   }
01085   DLIList<CubitFacet *> facet_list;
01086   CAST_LIST(facetentity_list, facet_list, CubitFacet);
01087 
01088   // set the control points on edges
01089 
01090   double edgelen;
01091   double compare_tol = 0.0;
01092   int index, ii;
01093   CubitFacetEdge *edge_ptr;
01094   for (ii=0; ii<numEdge; ii++)
01095   {
01096     edge_ptr = edge_array[ii];
01097     index = 3 * (ii * NUM_EDGE_CPTS);
01098     edge_ptr->set_control_points( &(edgeCtrlPts[index]) );
01099     if ( !specify_tol )
01100     {
01101         edgelen = edge_ptr->length();
01102         compare_tol += edgelen; 
01103     }
01104   }
01105   if ( specify_tol )
01106   {
01107       compare_tol = converge_tol;
01108   }
01109   else
01110   {
01111       compare_tol /= numEdge;
01112       compare_tol *= 1.0e-3;
01113   }
01114 
01115   // set the control points on triangles
01116 
01117   CubitFacet *facet_ptr;
01118   for (ii=0; ii<numTri; ii++)
01119   {
01120     facet_ptr = facet_array[ii];
01121     index = 3 * (ii * NUM_TRI_CPTS); 
01122     facet_ptr->set_control_points( &(faceCtrlPts[index]) );
01123   }
01124 
01125   // set the control points on quads
01126 
01127   CubitQuadFacet *qfacet_ptr;
01128   for (ii=0; ii<numQuad; ii++)
01129   {
01130     qfacet_ptr = qfacet_array[ii];
01131     index = 3 * (ii * NUM_QUAD_CPTS);
01132     qfacet_ptr->set_control_points( &(faceCtrlPts[index]) );
01133   }
01134 
01135   // Do the projections
01136 
01137   CubitFacet *last_facet = NULL;
01138   int interp_order = 4;
01139   CubitBoolean trim = CUBIT_TRUE;
01140   CubitBoolean outside;
01141   for (ii=0; ii<numLocs; ii++)
01142   {
01143     CubitVector cur_location(xyz[ii*3],
01144                              xyz[ii*3+1],
01145                              xyz[ii*3+2]);
01146     CubitVector eval_point;
01147     CubitVector *eval_point_ptr = &eval_point;
01148     CubitVector eval_normal;
01149     CubitVector *eval_normal_ptr = NULL;
01150     if (normal != NULL)
01151     {
01152       eval_normal_ptr = &eval_normal;
01153     }
01154     FacetEvalTool::project_to_facets(facet_list, last_facet,
01155                                      interp_order, compare_tol,
01156                                      cur_location, trim,
01157                                      &outside, eval_point_ptr,
01158                                      eval_normal_ptr);
01159     index = 3 * ii;
01160     xyzOnFace[index] = eval_point.x();
01161     xyzOnFace[index+1] = eval_point.y();
01162     xyzOnFace[index+2] = eval_point.z();
01163     if (normal != NULL)
01164     {
01165       normal[index] = eval_normal.x();
01166       normal[index+1] = eval_normal.y();
01167       normal[index+2] = eval_normal.z();
01168     }
01169   }
01170 
01171   // no derivatives yet!!
01172 
01173   
01174   // delete the temp arrays
01175 
01176   for (ii=0; ii<numQuad; ii++)
01177     delete qfacet_array[ii];
01178   for (ii=0; ii<numTri; ii++)
01179     delete facet_array[ii];
01180   for (ii=0; ii<numEdge; ii++)
01181     delete edge_array[ii];
01182   for (ii=0; ii<numVert; ii++)
01183     delete point_array[ii];
01184   delete [] point_array;
01185   delete [] edge_array;
01186   if (facet_array)
01187     delete [] facet_array; 
01188   if (qfacet_array)
01189     delete [] qfacet_array; 
01190 }
01191 
01192 //===========================================================================
01193 //  Function: projToBezierFaceFromNorms
01194 //  Purpose:  This is the same as the projToBezierFace except it differes by 
01195 //  the input arguments.  This function takes a list of normals and tangents 
01196 //  at the face vertices and computes bezier control points internally.  
01197 //  Normals and tangents should have been computed previously in 
01198 //  constructTriNormals or constructQuadNormals.  This function is not as 
01199 //  computationally efficient as evalBezierFace since it requires Bezier 
01200 //  control points to be computed as part of the call - however, it is more 
01201 //  memory efficient, requiring fewer variables to be stored with the calling 
01202 //  application. The following argument list describes only those arguments 
01203 //  that differ from projToBezierFace above.
01204 //
01205 //  Date:     12/08/2003
01206 //  Author:   sjowen
01207 //===========================================================================
01208 void projToBezierFaceFromNorms( int numFace, int numEdge, int numVert, 
01209                               int numVertPerFace, int* faceEdge, int* edgeVert, 
01210                               double* vert, double* vertNorms, double* edgeVertTang,
01211                               int numLocs, double* xyz, 
01212                               int specify_tol, double converge_tol,
01213                               double* xyzOnFace, double* normal, double* deriv )
01214 {
01215       // create arrays of facet entities from the arrays
01216 
01217   int numTri = 0;
01218   int numQuad = 0;
01219   CubitPoint **point_array = new CubitPoint * [numVert];
01220   CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
01221   CubitFacet **facet_array = NULL;
01222   int *triEdge = (numVertPerFace == 3) ? faceEdge : NULL;
01223   int *quadEdge = (numVertPerFace == 4) ? faceEdge : NULL;
01224   if (numVertPerFace == 3)
01225   {
01226     facet_array = new CubitFacet * [numFace];
01227     numTri = numFace;
01228   }
01229   CubitQuadFacet **qfacet_array = NULL;
01230   if (numVertPerFace == 4)
01231   {
01232     qfacet_array = new CubitQuadFacet * [numFace];
01233     numQuad = numFace;
01234   }
01235   DLIList<FacetEntity *> facetentity_list;
01236   CubitStatus stat = CUBIT_SUCCESS;
01237   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
01238                       edgeVert, vert, point_array, edge_array, facet_array,
01239                       qfacet_array, facetentity_list);
01240   if (CUBIT_SUCCESS != stat) {
01241     PRINT_ERROR("Failed to build facets.\n");
01242     return;
01243   }
01244   DLIList<CubitFacet *> facet_list;
01245   CAST_LIST(facetentity_list, facet_list, CubitFacet);
01246 
01247   // set the normals on triangle vertices
01248 
01249   int index, ii, jj;
01250   CubitPoint *pt;
01251   CubitFacet *facet_ptr;
01252   CubitVector pt_normal;
01253   for (ii=0; ii<numTri; ii++)
01254   {
01255     facet_ptr = facet_array[ii];
01256     for(jj=0; jj<3; jj++)
01257     {
01258       index = ii * 9 + (3 * jj);
01259       pt_normal.x( vertNorms[ index ] );
01260       pt_normal.y( vertNorms[ index + 1 ] );
01261       pt_normal.z( vertNorms[ index + 2 ] );
01262       pt = facet_ptr->point(jj);
01263       TDFacetBoundaryPoint::add_facet_boundary_point( pt, facet_ptr, pt_normal );
01264     }
01265   }
01266 
01267   // set the normals on quad vertices
01268 
01269   CubitQuadFacet *qfacet_ptr;
01270   for (ii=0; ii<numQuad; ii++)
01271   {
01272     qfacet_ptr = qfacet_array[ii];
01273     for(jj=0; jj<4; jj++)
01274     {
01275       index = ii * 12 + (3 * jj);
01276       pt_normal.x( vertNorms[ index ] );
01277       pt_normal.y( vertNorms[ index + 1 ] );
01278       pt_normal.z( vertNorms[ index + 2 ] );
01279       pt = qfacet_ptr->point(jj);
01280       TDFacetBoundaryPoint::add_facet_boundary_point( pt, qfacet_ptr, pt_normal );
01281     }
01282   }
01283 
01284   // set the control points on edges
01285 
01286   double edgelen;
01287   double compare_tol = 0.0;
01288   CubitFacetEdge *edge_ptr;
01289   CubitVector N0, N1, T0, T1, P0, P1;
01290   CubitVector Pi[3];
01291   CubitPoint *pt0, *pt1;
01292   for (ii=0; ii<numEdge; ii++)
01293   {
01294     index = 6 * ii;
01295     edge_ptr = edge_array[ii];
01296     pt0 = edge_ptr->point(0);
01297     pt1 = edge_ptr->point(1);
01298     P0 = pt0->coordinates();
01299     P1 = pt1->coordinates();
01300     N0 = pt0->normal( edge_ptr );
01301     N1 = pt1->normal( edge_ptr );
01302     T0.x( edgeVertTang[ index ] );
01303     T0.y( edgeVertTang[ index + 1] );
01304     T0.z( edgeVertTang[ index + 2] );
01305     T1.x( edgeVertTang[ index + 3] );
01306     T1.y( edgeVertTang[ index + 4] );
01307     T1.z( edgeVertTang[ index + 5] );
01308     FacetEvalTool::init_edge_control_points(P0, P1, N0, N1, T0, T1, Pi);
01309     edge_ptr->control_points( Pi, 4 );
01310     if ( !specify_tol )
01311     {
01312         edgelen = edge_ptr->length();
01313         compare_tol += edgelen; 
01314     }
01315   }
01316   if ( specify_tol )
01317   {
01318       compare_tol = converge_tol;
01319   }
01320   else
01321   {
01322       compare_tol /= numEdge;
01323       compare_tol *= 1.0e-3;
01324   }
01325 
01326   // set the control points on triangles
01327 
01328   for (ii=0; ii<numTri; ii++)
01329   {
01330     facet_ptr = facet_array[ii];
01331     facet_ptr->init_patch();
01332   }
01333 
01334     // set the control points on quads
01335 
01336   for (ii=0; ii<numQuad; ii++)
01337   {
01338     qfacet_ptr = qfacet_array[ii];
01339     qfacet_ptr->init_patch();
01340   }
01341 
01342     // Do the projections
01343 
01344   CubitFacet *last_facet = NULL;
01345   int interp_order = 4;
01346   CubitBoolean trim = CUBIT_TRUE;
01347   CubitBoolean outside;
01348   
01349   for (ii=0; ii<numLocs; ii++)
01350   {
01351     CubitVector cur_location(xyz[ii*3],
01352                              xyz[ii*3+1],
01353                              xyz[ii*3+2]);
01354     CubitVector eval_point;
01355     CubitVector *eval_point_ptr = &eval_point;
01356     CubitVector eval_normal;
01357     CubitVector *eval_normal_ptr = NULL;
01358     if (normal != NULL)
01359     {
01360       eval_normal_ptr = &eval_normal;
01361     }
01362     FacetEvalTool::project_to_facets(facet_list, last_facet,
01363                                      interp_order, compare_tol,
01364                                      cur_location, trim,
01365                                      &outside, eval_point_ptr,
01366                                      eval_normal_ptr);
01367       //double dist = cur_location.distance_between(eval_point);
01368     index = 3 * ii;
01369     xyzOnFace[index] = eval_point.x();
01370     xyzOnFace[index+1] = eval_point.y();
01371     xyzOnFace[index+2] = eval_point.z();
01372     if (normal != NULL)
01373     {
01374       normal[index] = eval_normal.x();
01375       normal[index+1] = eval_normal.y();
01376       normal[index+2] = eval_normal.z();
01377     }
01378   }
01379 
01380   // no derivatives yet!!
01381 
01382   // delete the temp arrays
01383 
01384   for (ii=0; ii<numQuad; ii++)
01385     delete qfacet_array[ii];
01386   for (ii=0; ii<numTri; ii++)
01387     delete facet_array[ii];
01388   for (ii=0; ii<numEdge; ii++)
01389     delete edge_array[ii];
01390   for (ii=0; ii<numVert; ii++)
01391     delete point_array[ii];
01392   delete [] point_array;
01393   delete [] edge_array;
01394   if (facet_array)
01395     delete [] facet_array; 
01396   if (qfacet_array)
01397     delete [] qfacet_array; 
01398 }
01399 
01400 //---------------------------------------------------------------------------
01401 // Functions for debugging
01402 //---------------------------------------------------------------------------
01403 
01404 //===========================================================================
01405 //  Function: dumpMesh
01406 //  Purpose:  dump the face mesh to a file that can be read by CUBIT
01407 //            Use the same definition of parameters as used with 
01408 //            resolveFaceVectors
01409 //  Date:     10/28/2002
01410 //  Author:   sjowen
01411 //===========================================================================
01412 void dumpMesh(const char *fileName, int includeResults, double angle,
01413               int numTri, int numQuad, int numEdge,
01414               int numVert, int* triEdge, int* quadEdge,
01415               int* edgeVert, double* vert,
01416               double* edgeCtrlPts, double* triCtrlPts,
01417               double* quadCtrlPts)
01418 {
01419   FILE *fp = fopen(fileName, "w");
01420   assert(fp != NULL); // couldn't open file for writing
01421   
01422   // write the header info
01423 
01424   fprintf(fp, "Cholla version 1.0\n");
01425   time_stamp( fp );
01426   fprintf(fp, "%d %d %d %d %d\n", numTri, numQuad, numEdge, numVert, includeResults);
01427   fprintf(fp, "%f\n", angle);
01428   
01429   // faceEdgeBegin
01430 
01431   int ii;
01432 
01433   // faceEdge
01434 
01435   fprintf(fp, "triEdge\n");
01436   for(ii=0; ii<numTri; ii++)
01437   {
01438     fprintf(fp, "%d %d %d\n", triEdge[ii*3], triEdge[ii*3+1], triEdge[ii*3+2]);
01439   }
01440 
01441   fprintf(fp, "quadEdge\n");
01442   for(ii=0; ii<numQuad; ii++)
01443   {
01444     fprintf(fp, "%d %d %d %d\n", quadEdge[ii*4], quadEdge[ii*4+1],
01445                                  quadEdge[ii*4+2], quadEdge[ii*4+3]);
01446   }
01447 
01448   // edgeVert
01449 
01450   fprintf(fp, "edgeVert\n");
01451   for (ii=0; ii<numEdge; ii++)
01452   {
01453     fprintf(fp, "%d %d\n", edgeVert[ii*2], edgeVert[ii*2+1]);
01454   }
01455 
01456   // vert
01457 
01458   fprintf(fp, "vert\n");
01459   for (ii=0; ii<numVert; ii++)
01460   {
01461     fprintf(fp, "%f %f %f\n", vert[ii*3], vert[ii*3+1], vert[ii*3+2]);
01462   }
01463 
01464   if (includeResults != 0)
01465   {
01466      // edge control points
01467 
01468     int jj;
01469     fprintf(fp, "edgeCtrlPts\n");
01470     for(ii=0; ii<numEdge; ii++)
01471     {
01472       for (jj=0; jj<NUM_EDGE_CPTS; jj++)
01473       {
01474         fprintf(fp, "%f %f %f\n", edgeCtrlPts[ii*jj*3], 
01475                                      edgeCtrlPts[ii*jj*3+1], 
01476                                      edgeCtrlPts[ii*jj*3+2]);
01477       }
01478     }
01479 
01480     // triangle control points
01481 
01482     fprintf(fp, "triCtrlPts\n");
01483     for(ii=0; ii<numTri; ii++)
01484     {
01485       for (jj=0; jj<NUM_TRI_CPTS; jj++)
01486       {
01487         fprintf(fp, "%f %f %f\n", triCtrlPts[ii*jj*3], 
01488                                      triCtrlPts[ii*jj*3+1], 
01489                                      triCtrlPts[ii*jj*3+2]);
01490       }
01491     }
01492 
01493     // quad control points
01494 
01495     fprintf(fp, "quadCtrlPts\n");
01496     for(ii=0; ii<numQuad; ii++)
01497     {
01498       for (jj=0; jj<NUM_QUAD_CPTS; jj++)
01499       {
01500         fprintf(fp, "%f %f %f\n", quadCtrlPts[ii*jj*3], 
01501                                      quadCtrlPts[ii*jj*3+1], 
01502                                      quadCtrlPts[ii*jj*3+2]);
01503       }
01504     }
01505   }
01506   fclose(fp);
01507 }
01508 
01509 //===========================================================================
01510 //  Function: dumpResults
01511 //  Purpose:  dump only the results to a file
01512 //  Date:     10/28/2002
01513 //  Author:   sjowen
01514 //===========================================================================
01515 void dumpResults(const char *fileName, int numEdge, int numTri, int numQuad, 
01516                  double* edgeCtrlPts, double* triCtrlPts, 
01517                  double* quadCtrlPts )
01518 {
01519   FILE *fp = fopen( fileName, "w" );
01520   assert(fp != NULL);
01521 
01522     // write the header info
01523 
01524   fprintf(fp, "Cholla version 1.0 Results\n");
01525   time_stamp( fp );
01526   fprintf(fp, "%d %d %d\n", numEdge, numTri, numQuad );
01527 
01528   int ii,jj;
01529   fprintf(fp, "edgeCtrlPts\n");
01530   for(ii=0; ii<numEdge; ii++)
01531   {
01532     for (jj=0; jj<NUM_EDGE_CPTS; jj++)
01533     {
01534       fprintf(fp, "%f %f %f\n", edgeCtrlPts[ii*jj*3], 
01535                                    edgeCtrlPts[ii*jj*3+1], 
01536                                    edgeCtrlPts[ii*jj*3+2]);
01537     }
01538   }
01539 
01540   // triangle control points
01541 
01542   fprintf(fp, "triCtrlPts\n");
01543   for(ii=0; ii<numTri; ii++)
01544   {
01545     for (jj=0; jj<NUM_TRI_CPTS; jj++)
01546     {
01547       fprintf(fp, "%f %f %f\n", triCtrlPts[ii*jj*3], 
01548                                    triCtrlPts[ii*jj*3+1], 
01549                                    triCtrlPts[ii*jj*3+2]);
01550     }
01551   }
01552 
01553   // quad control points
01554 
01555   fprintf(fp, "quadCtrlPts\n");
01556   for(ii=0; ii<numQuad; ii++)
01557   {
01558     for (jj=0; jj<NUM_QUAD_CPTS; jj++)
01559     {
01560       fprintf(fp, "%f %f %f\n", quadCtrlPts[ii*jj*3], 
01561                                    quadCtrlPts[ii*jj*3+1], 
01562                                    quadCtrlPts[ii*jj*3+2]);
01563     }
01564   }
01565 
01566   fclose(fp);
01567 }
01568 
01569 //===========================================================================
01570 //  Function: importMesh
01571 //  Purpose:  import the face mesh from the specified file
01572 //            allocate the arrays and pass them back
01573 //  Date:     11/28/2002
01574 //  Author:   sjowen
01575 //===========================================================================
01576 void importMesh(const char *fileName, int *includeResults,
01577                 double *angle, int *numTri, int *numQuad, int *numEdge,
01578                 int *numVert, int** triEdge, int** quadEdge,
01579                 int** edgeVert, double** vert,
01580                 double** edgeCtrlPts, double** triCtrlPts,
01581                 double** quadCtrlPts)
01582 {
01583   FILE *fp = fopen(fileName, "r");
01584   assert(fp != NULL);  // couldn't open file for reading
01585   
01586   // read the header info
01587 
01588   char version[128];
01589   char* s = fgets(version, 128, fp);
01590   if (NULL == s) {
01591     PRINT_ERROR("Format error in file %s\n", fileName);
01592     fclose(fp);
01593     return;
01594   }
01595   assert(strcmp(version, "Cholla version 1.0\n") == 0);  // don't recognize file
01596   char filetime[128];
01597   s = fgets(filetime, 128, fp); 
01598   if (NULL == s) {
01599     PRINT_ERROR("Format error in file %s\n", fileName);
01600     fclose(fp);
01601     return;
01602   }
01603   
01604   char line[256];
01605   int num_tri;
01606   int num_quad;
01607   int num_edge;
01608   int num_vert;
01609   s = fgets(line,256,fp);
01610   if (NULL == s) {
01611     PRINT_ERROR("Format error in file %s\n", fileName);
01612     fclose(fp);
01613     return;
01614   }
01615   sscanf(line, "%d %d %d %d %d\n", &num_tri, &num_quad, &num_edge, &num_vert, includeResults); 
01616   *numTri = num_tri;
01617   *numQuad = num_quad;
01618   *numEdge = num_edge;
01619   *numVert = num_vert;
01620   int nread = fscanf(fp, "%lf\n", angle);
01621   if (1 != nread) {
01622     PRINT_ERROR("Format error in file %s\n", fileName);
01623     fclose(fp);
01624     return;
01625   }
01626 
01627   // triEdge
01628 
01629   *triEdge = new int [num_tri*3];
01630   int *tri_edge = *triEdge;
01631   char array_name[128];
01632   nread = fscanf(fp, "%s\n", array_name);
01633   if (1 != nread) {
01634     PRINT_ERROR("Format error in file %s\n", fileName);
01635     fclose(fp);
01636     return;
01637   }
01638   assert(strcmp(array_name, "triEdge") == 0);  // check start of array
01639   int ii;
01640   for(ii=0; ii<num_tri; ii++)
01641   {
01642     nread = fscanf(fp, "%d %d %d\n", &tri_edge[ii*3], &tri_edge[ii*3+1], &tri_edge[ii*3+2]);
01643     if (3 != nread) {
01644       PRINT_ERROR("Format error in file %s\n", fileName);
01645       fclose(fp);
01646       return;
01647     }
01648   }
01649 
01650   // quadEdge
01651 
01652   *quadEdge = new int [num_quad*4];
01653   int *quad_edge = *quadEdge;
01654   nread = fscanf(fp, "%s\n", array_name);
01655   if (1 != nread) {
01656     PRINT_ERROR("Format error in file %s\n", fileName);
01657     fclose(fp);
01658     return;
01659   }
01660   assert(strcmp(array_name, "quadEdge") == 0);  // check start of array
01661   for(ii=0; ii<num_quad; ii++)
01662   {
01663     nread = fscanf(fp, "%d %d %d %d\n", &quad_edge[ii*4], &quad_edge[ii*4+1], 
01664                                 &quad_edge[ii*4+2], &quad_edge[ii*4+3]);
01665     if (4 != nread) {
01666       PRINT_ERROR("Format error in file %s\n", fileName);
01667       fclose(fp);
01668       return;
01669     }
01670   }
01671 
01672   // edgeVert
01673 
01674   *edgeVert = new int [2*num_edge];
01675   int *edge_vert = *edgeVert;
01676   nread = fscanf(fp, "%s\n", array_name);
01677   if (1 != nread) {
01678     PRINT_ERROR("Format error in file %s\n", fileName);
01679     fclose(fp);
01680     return;
01681   }
01682   assert(strcmp(array_name, "edgeVert") == 0);  // check start of array
01683   for (ii=0; ii<num_edge; ii++)
01684   {
01685     nread = fscanf(fp, "%d %d\n", &edge_vert[ii*2], &edge_vert[ii*2+1]);
01686     if (2 != nread) {
01687       PRINT_ERROR("Format error in file %s\n", fileName);
01688       fclose(fp);
01689       return;
01690     }
01691   }
01692 
01693   // vert
01694 
01695   *vert = new double [3*num_vert];
01696   double *my_vert = *vert;
01697   nread = fscanf(fp, "%s\n", array_name);
01698   if (1 != nread) {
01699     PRINT_ERROR("Format error in file %s\n", fileName);
01700     fclose(fp);
01701     return;
01702   }
01703   assert(strcmp(array_name, "vert") == 0);  // check start of array
01704   for (ii=0; ii<num_vert; ii++)
01705   {
01706     nread = fscanf(fp, "%lf %lf %lf\n", &my_vert[ii*3], &my_vert[ii*3+1], &my_vert[ii*3+2]);
01707     if (3 != nread) {
01708       PRINT_ERROR("Format error in file %s\n", fileName);
01709       fclose(fp);
01710       return;
01711     }
01712   }
01713 
01714   *edgeCtrlPts = new double [3*num_edge*NUM_EDGE_CPTS];
01715   if (numTri != 0)
01716     *triCtrlPts = new double [3*num_tri*NUM_TRI_CPTS];
01717   if (numQuad != 0)
01718     *quadCtrlPts = new double [3*num_quad*NUM_QUAD_CPTS];
01719 
01720   double *edge_ctrl_pts = *edgeCtrlPts;
01721   double *tri_ctrl_pts = *triCtrlPts;
01722   double *quad_ctrl_pts = *quadCtrlPts;
01723   if (*includeResults != 0)
01724   {
01725      // edge control points
01726     
01727     nread = fscanf(fp, "%s\n", array_name);
01728     if (1 != nread) {
01729       PRINT_ERROR("Format error in file %s\n", fileName);
01730       fclose(fp);
01731       return;
01732     }
01733     assert(strcmp(array_name, "edgeCtrlPts") == 0);  // check start of array
01734     for(ii=0; ii<num_edge*NUM_EDGE_CPTS; ii++)
01735     {
01736       nread = fscanf(fp, "%lf %lf %lf\n", &edge_ctrl_pts[ii*3], 
01737                                   &edge_ctrl_pts[ii*3+1], 
01738                                   &edge_ctrl_pts[ii*3+2]);
01739       if (3 != nread) {
01740         PRINT_ERROR("Format error in file %s\n", fileName);
01741         fclose(fp);
01742         return;
01743       }
01744     }
01745 
01746      // triangle control points
01747     
01748     nread = fscanf(fp, "%s\n", array_name);
01749     if (1 != nread) {
01750       PRINT_ERROR("Format error in file %s\n", fileName);
01751       fclose(fp);
01752       return;
01753     }
01754     assert(strcmp(array_name, "triCtrlPts") == 0);  // check start of array
01755     for(ii=0; ii<num_tri*NUM_TRI_CPTS; ii++)
01756     {
01757       nread = fscanf(fp, "%lf %lf %lf\n", &tri_ctrl_pts[ii*3], 
01758                                   &tri_ctrl_pts[ii*3+1], 
01759                                   &tri_ctrl_pts[ii*3+2]);
01760       if (3 != nread) {
01761         PRINT_ERROR("Format error in file %s\n", fileName);
01762         fclose(fp);
01763         return;
01764       }
01765     }
01766 
01767     // quad control points
01768     
01769     nread = fscanf(fp, "%s\n", array_name);
01770     if (1 != nread) {
01771       PRINT_ERROR("Format error in file %s\n", fileName);
01772       fclose(fp);
01773       return;
01774     }
01775     assert(strcmp(array_name, "quadCtrlPts") == 0);  // check start of array
01776     for(ii=0; ii<num_quad*NUM_QUAD_CPTS; ii++)
01777     {
01778       nread = fscanf(fp, "%lf %lf %lf\n", &quad_ctrl_pts[ii*3], 
01779                                   &quad_ctrl_pts[ii*3+1], 
01780                                   &quad_ctrl_pts[ii*3+2]);
01781       if (3 != nread) {
01782         PRINT_ERROR("Format error in file %s\n", fileName);
01783         fclose(fp);
01784         return;
01785       }
01786     }
01787 
01788   }
01789   fclose(fp);
01790 }
01791 
01792 //===========================================================================
01793 //  Function: importResults
01794 //  Purpose:  import only the results to from a file
01795 //  Note:     allocates the normal and tangent arrays
01796 //  Date:     10/28/2002
01797 //  Author:   sjowen
01798 //===========================================================================
01799 void importResults(const char *fileName, int *numEdge, int *numTri, 
01800                    int *numQuad, double** edgeCtrlPts, double** triCtrlPts, 
01801                    double** quadCtrlPts )
01802 {
01803   FILE *fp = fopen(fileName, "r");
01804   assert(fp != NULL);  // couldn't open file for reading
01805   
01806   // read the header info
01807 
01808   char version[128];
01809   char* s = fgets(version, 128, fp);
01810   if (NULL == s) {
01811     PRINT_ERROR("Format error in file %s\n", fileName);
01812     fclose(fp);
01813     return;
01814   }
01815   assert(strcmp(version, "Cholla version 1.0 Results\n") == 0);  // don't recognize file
01816   char filetime[128];
01817   s = fgets(filetime, 128, fp);
01818   if (NULL == s) {
01819     PRINT_ERROR("Format error in mesh file %s\n", fileName);
01820     fclose(fp);
01821     return;
01822   }
01823   int nread = fscanf(fp, "%d %d %d\n", numEdge, numTri, numQuad);
01824   if (3 != nread) {
01825     PRINT_ERROR("Format error in file %s\n", fileName);
01826     fclose(fp);
01827     return;
01828   }
01829   int num_edge = *numEdge;
01830   int num_tri = *numTri;
01831   int num_quad = *numQuad;
01832   
01833   *edgeCtrlPts = new double [3*num_edge*NUM_EDGE_CPTS];
01834   if (numTri != 0)
01835     *triCtrlPts = new double [3*num_tri*NUM_TRI_CPTS];
01836   if (numQuad != 0)
01837     *quadCtrlPts = new double [3*num_quad*NUM_QUAD_CPTS];
01838 
01839   double *edge_ctrl_pts = *edgeCtrlPts;
01840   double *tri_ctrl_pts = *triCtrlPts;
01841   double *quad_ctrl_pts = *quadCtrlPts;
01842 
01843    // edge control points
01844   
01845   int ii;
01846   char array_name[128];
01847   nread = fscanf(fp, "%s\n", array_name);
01848   if (1 != nread) {
01849     PRINT_ERROR("Format error in file %s\n", fileName);
01850     fclose(fp);
01851     return;
01852   }
01853   assert(strcmp(array_name, "edgeCtrlPts") == 0);  // check start of array
01854   for(ii=0; ii<num_edge*NUM_EDGE_CPTS; ii++)
01855   {
01856     nread = fscanf(fp, "%lf %lf %lf\n", &edge_ctrl_pts[ii*3], 
01857                                 &edge_ctrl_pts[ii*3+1], 
01858                                 &edge_ctrl_pts[ii*3+2]);
01859     if (3 != nread) {
01860       PRINT_ERROR("Format error in file %s\n", fileName);
01861       fclose(fp);
01862       return;
01863     }
01864   }
01865 
01866    // triangle control points
01867   
01868   nread = fscanf(fp, "%s\n", array_name);
01869   if (1 != nread) {
01870     PRINT_ERROR("Format error in file %s\n", fileName);
01871     fclose(fp);
01872     return;
01873   }
01874   assert(strcmp(array_name, "triCtrlPts") == 0);  // check start of array
01875   for(ii=0; ii<num_tri*NUM_TRI_CPTS; ii++)
01876   {
01877     nread = fscanf(fp, "%lf %lf %lf\n", &tri_ctrl_pts[ii*3], 
01878                                 &tri_ctrl_pts[ii*3+1], 
01879                                 &tri_ctrl_pts[ii*3+2]);
01880     if (3 != nread) {
01881       PRINT_ERROR("Format error in file %s\n", fileName);
01882       fclose(fp);
01883       return;
01884     }
01885   }
01886 
01887   // quad control points
01888   
01889   nread = fscanf(fp, "%s\n", array_name);
01890   if (1 != nread) {
01891     PRINT_ERROR("Format error in file %s\n", fileName);
01892     fclose(fp);
01893     return;
01894   }
01895   assert(strcmp(array_name, "quadCtrlPts") == 0);  // check start of array
01896   for(ii=0; ii<num_tri*NUM_QUAD_CPTS; ii++)
01897   {
01898     nread = fscanf(fp, "%lf %lf %lf\n", &quad_ctrl_pts[ii*3], 
01899                                 &quad_ctrl_pts[ii*3+1], 
01900                                 &quad_ctrl_pts[ii*3+2]);
01901     if (3 != nread) {
01902       PRINT_ERROR("Format error in file %s\n", fileName);
01903       fclose(fp);
01904       return;
01905     }
01906   }
01907 
01908   fclose(fp);
01909 
01910 }
01911 
01912 //===========================================================================
01913 //  Function: time_stamp
01914 //  Purpose:  write the current time to a file
01915 //  Date:     11/28/2002
01916 //  Author:   sjowen
01917 //===========================================================================
01918 static void time_stamp( FILE *fp )
01919 {
01920   struct tm *newtime;
01921   bool am = true;
01922   time_t long_time;
01923 
01924   time( &long_time );                /* Get time as long integer. */
01925   newtime = localtime( &long_time ); /* Convert to local time. */
01926 
01927   if( newtime->tm_hour > 12 )        /* Set up extension. */
01928     am = false;
01929   if( newtime->tm_hour > 12 )        /* Convert from 24-hour */
01930           newtime->tm_hour -= 12;    /*   to 12-hour clock.  */
01931   if( newtime->tm_hour == 0 )        /*Set hour to 12 if midnight. */
01932           newtime->tm_hour = 12;
01933 
01934   fprintf( fp, "%.19s %s\n", asctime( newtime ), am ? "AM" : "PM" );
01935 }
01936 
01937 //===========================================================================
01938 //  Function: evalBezierEdge
01939 //  Purpose:  Evaluate a quartic Bezier curve at a specified parametric
01940 //  location given its control points. Evaluates location, and/or tangent.
01941 
01942 //  Date:     01/20/2004
01943 //  Author:   jfowler
01944 //===========================================================================
01945 void evalBezierEdge( int numEdge, int numVert, int* edgeVert,
01946                      double* vert, double* edgeCtrlPts, 
01947                      int numLocs, double* paramLocation, 
01948                      double* location, double* tangent )
01949 {
01950 CubitPoint **point_array = new CubitPoint * [numVert];
01951 CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
01952 CubitFacet **facet_array = NULL;
01953 CubitQuadFacet **qfacet_array = NULL;
01954 DLIList<FacetEntity *> facet_list;
01955 int ii, index, numTri, numQuad;
01956 int *triEdge, *quadEdge;
01957   triEdge = quadEdge = NULL;
01958   numTri = numQuad = 0;
01959   CubitStatus stat = CUBIT_SUCCESS;
01960   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
01961                       edgeVert, vert, point_array, edge_array, facet_array,
01962                       qfacet_array, facet_list);
01963   if (CUBIT_SUCCESS != stat) {
01964     PRINT_ERROR("Failed to build facets.\n");
01965     return;
01966   }
01967 
01968   CubitFacetEdge *edge_ptr;
01969   //  Set the control points on the edges.
01970   for (ii=0; ii<numEdge; ii++)
01971   {
01972     edge_ptr = edge_array[ii];
01973     index = 3 * (ii * NUM_EDGE_CPTS);
01974     edge_ptr->set_control_points( &(edgeCtrlPts[index]) );
01975   } 
01976 
01977   // do the evaluations
01978   int jj, tindex;
01979   index = tindex = 0;
01980     
01981     for ( ii = 0; ii < numEdge; ii++ ) {
01982     for ( jj = 0; jj < numLocs; jj++ ) {
01983       
01984       CubitVector eval_point;
01985       CubitVector *eval_point_ptr = &eval_point;
01986       CubitVector eval_tangent;
01987       CubitVector *eval_tangent_ptr = NULL;
01988       if ( tangent != NULL ) 
01989         eval_tangent_ptr = &eval_tangent;
01990       edge_ptr = edge_array[ii];
01991       edge_ptr->evaluate_single(paramLocation[jj],eval_point_ptr);      
01992       location[index++] = eval_point.x();
01993       location[index++] = eval_point.y();
01994       location[index++] = eval_point.z();
01995       if ( tangent != NULL ) {  
01996         edge_ptr->evaluate_single_tangent(paramLocation[jj],eval_tangent_ptr);      
01997         eval_tangent.normalize(); 
01998         tangent[tindex++] = eval_tangent.x();
01999         tangent[tindex++] = eval_tangent.y();
02000         tangent[tindex++] = eval_tangent.z();             
02001       }      
02002     }  
02003   }
02004   //  delete the temp arrays
02005 
02006   for (ii=0; ii<numEdge; ii++)
02007     delete edge_array[ii];
02008   for (ii=0; ii<numVert; ii++)
02009     delete point_array[ii];
02010   delete [] point_array;
02011   delete [] edge_array;
02012 
02013 }
02014 
02015 //===========================================================================
02016 //  Function: evalBezierEdgeFromTans
02017 //  Purpose:  Evaluate a quartic Bezier curve at a specified parametric
02018 //  location given its tangents at the end-points. Evaluates location, 
02019 //  and/or tangent.
02020 //  Date:     01/20/2004
02021 //  Author:   jfowler
02022 //===========================================================================
02023 void evalBezierEdgeFromTans( int numEdge, int numVert, int* edgeVert, 
02024                               double* vert, double* edgeVertTang,
02025                               int numLocs, double* paramLocation, 
02026                               double* location, double* tangent )
02027 {
02028 CubitPoint **point_array = new CubitPoint * [numVert];
02029 CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
02030 CubitFacet **facet_array = NULL;
02031 CubitQuadFacet **qfacet_array = NULL;
02032 DLIList<FacetEntity *> facet_list;
02033 int ii, index, numTri, numQuad;
02034 int *triEdge, *quadEdge;
02035   triEdge = quadEdge = NULL;
02036   numTri = numQuad = 0;
02037   CubitStatus stat = CUBIT_SUCCESS;
02038   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
02039                       edgeVert, vert, point_array, edge_array, facet_array,
02040                       qfacet_array, facet_list);
02041   if (CUBIT_SUCCESS != stat) {
02042     PRINT_ERROR("Failed to build facets.\n");
02043     return;
02044   }
02045 
02046   // set the control points on edges
02047 
02048   int mydebug = 0;
02049   FILE *fp = NULL;
02050   if (mydebug)
02051     fp = fopen("edges.txt", "w");
02052   CubitFacetEdge *edge_ptr;
02053   CubitVector T0, T1, P0, P1;
02054   CubitVector Pi[3];
02055   CubitPoint *pt0, *pt1;
02056   for (ii=0; ii<numEdge; ii++)
02057   {
02058     index = 6 * ii;
02059     edge_ptr = edge_array[ii];
02060     pt0 = edge_ptr->point(0);
02061     pt1 = edge_ptr->point(1);
02062     P0 = pt0->coordinates();
02063     P1 = pt1->coordinates();
02064     T0.x( edgeVertTang[ index ] );
02065     T0.y( edgeVertTang[ index + 1] );
02066     T0.z( edgeVertTang[ index + 2] );
02067     T1.x( edgeVertTang[ index + 3] );
02068     T1.y( edgeVertTang[ index + 4] );
02069     T1.z( edgeVertTang[ index + 5] );
02070     FacetEvalTool::init_edge_control_points_single(P0, P1, T0, T1, Pi);
02071     if (mydebug)
02072     {
02073       int kk;
02074       for(kk=0; kk<3; kk++);
02075       fprintf(fp, "(%d) %.6f %.6f %.6f\n", edge_ptr->id(), P0.x(), P0.y(), P0.z());
02076       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[0].x(), Pi[0].y(), Pi[0].z());
02077       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[1].x(), Pi[1].y(), Pi[1].z());
02078       fprintf(fp, "    %.6f %.6f %.6f\n", Pi[2].x(), Pi[2].y(), Pi[2].z());
02079       fprintf(fp, "    %.6f %.6f %.6f\n", P1.x(), P1.y(), P1.z());
02080     }
02081     edge_ptr->control_points( Pi, 4 );
02082   } 
02083   if (mydebug)
02084     fclose(fp);
02085 
02086   // do the evaluations
02087   int jj, tindex;
02088   index = tindex = 0;
02089     
02090   for ( ii = 0; ii < numEdge; ii++ ) {
02091     for ( jj = 0; jj < numLocs; jj++ ) {
02092       
02093       CubitVector eval_point;
02094       CubitVector *eval_point_ptr = &eval_point;
02095       CubitVector eval_tangent;
02096       CubitVector *eval_tangent_ptr = NULL;
02097       if ( tangent != NULL ) 
02098         eval_tangent_ptr = &eval_tangent;
02099       edge_ptr = edge_array[ii];
02100       edge_ptr->evaluate_single(paramLocation[jj],eval_point_ptr);      
02101       location[index++] = eval_point.x();
02102       location[index++] = eval_point.y();
02103       location[index++] = eval_point.z();
02104       if ( tangent != NULL ) {  
02105         edge_ptr->evaluate_single_tangent(paramLocation[jj],eval_tangent_ptr);      
02106         eval_tangent.normalize(); 
02107         tangent[tindex++] = eval_tangent.x();
02108         tangent[tindex++] = eval_tangent.y();
02109         tangent[tindex++] = eval_tangent.z();             
02110       }      
02111     }  
02112   }
02113   //  delete the temp arrays
02114 
02115   for (ii=0; ii<numEdge; ii++)
02116     delete edge_array[ii];
02117   for (ii=0; ii<numVert; ii++)
02118     delete point_array[ii];
02119   delete [] point_array;
02120   delete [] edge_array;
02121 }
02122 
02123 //===========================================================================
02124 //  Function: projBezierEdgeFromTans
02125 //  Purpose:  Project a point to an edge given tangents at the endpoints.
02126 //            The routine only finds one point.  Rarely there might be
02127 //            two or more closest points.
02128 //  Date:     01/20/2004             
02129 //  Author:   jfowler
02130 //===========================================================================
02131 void projToBezierEdgeFromTans( int numEdge, int numVert, int* edgeVert, 
02132                               double* vert, double* edgeVertTang,
02133                               int numLocs, double* xyz, 
02134                               double* xyzOnEdge, double* tangent )
02135 {
02136 CubitPoint **point_array = new CubitPoint * [numVert];
02137 CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
02138 CubitFacet **facet_array = NULL;
02139 CubitQuadFacet **qfacet_array = NULL;
02140 DLIList<FacetEntity *> facet_list;
02141 int ii, index, numTri, numQuad;
02142 int *triEdge, *quadEdge;
02143   triEdge = quadEdge = NULL;
02144   numTri = numQuad = 0;
02145   CubitStatus stat = CUBIT_SUCCESS;
02146   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
02147                       edgeVert, vert, point_array, edge_array, facet_array,
02148                       qfacet_array, facet_list);
02149   if (CUBIT_SUCCESS != stat) {
02150     PRINT_ERROR("Failed to build facets.\n");
02151     return;
02152   }
02153 
02154   // set the control points on edges
02155 
02156   CubitFacetEdge *edge_ptr;
02157   CubitVector T0, T1, P0, P1;
02158   CubitVector Pi[3];
02159   CubitPoint *pt0, *pt1;
02160   for (ii=0; ii<numEdge; ii++)
02161   {
02162     index = 6 * ii;
02163     edge_ptr = edge_array[ii];
02164     pt0 = edge_ptr->point(0);
02165     pt1 = edge_ptr->point(1);
02166     P0 = pt0->coordinates();
02167     P1 = pt1->coordinates();
02168     T0.x( edgeVertTang[ index ] );
02169     T0.y( edgeVertTang[ index + 1] );
02170     T0.z( edgeVertTang[ index + 2] );
02171     T1.x( edgeVertTang[ index + 3] );
02172     T1.y( edgeVertTang[ index + 4] );
02173     T1.z( edgeVertTang[ index + 5] );
02174     FacetEvalTool::init_edge_control_points_single(P0, P1, T0, T1, Pi);
02175 
02176     edge_ptr->control_points( Pi, 4 );
02177   } 
02178 
02179   // do the projections
02180   int jj, kk, xyzindex, tindex;
02181   index = tindex = 0;
02182   double dsquared, dsquared_test, dderiv, dderiv2, t, t_min = 0.0;
02183   CubitVector xyz_pt;//, second_d;
02184   t = 0.0;
02185   dsquared_test =  CUBIT_DBL_MAX;  
02186   for ( ii = 0; ii < numEdge; ii++ ) {
02187     edge_ptr = edge_array[ii];
02188     for ( jj = 0; jj < numLocs; jj++ ) {
02189       xyzindex = 3*jj;
02190       xyz_pt.x(xyz[xyzindex]);
02191       xyz_pt.y(xyz[xyzindex+1]);
02192       xyz_pt.z(xyz[xyzindex+2]);    
02193       t = -1.0;
02194       dsquared_test =  CUBIT_DBL_MAX;  
02195       //  Get square of distance from xyz to curve at intervals of 0.2
02196       //  in the parameter t.  Save the closest distance found as a 
02197       //  starting point for Newton-Raphson iterations.  It was found 
02198       //  that an interval of 0.5 was too coarse and that the N-R
02199       //  sometimes converged on a false value.
02200       CubitVector eval_point;
02201       CubitVector *eval_point_ptr = &eval_point;
02202       CubitVector eval_tangent;
02203       CubitVector *eval_tangent_ptr = &eval_tangent;
02204       CubitVector second_d;
02205       CubitVector *second_d_ptr = &second_d;
02206       for ( kk = 0; kk < 11; kk++ ) {
02207         edge_ptr->evaluate_single(t,eval_point_ptr);      
02208         dsquared = ( (xyz_pt.x()-eval_point.x())*(xyz_pt.x()-eval_point.x()) +
02209                      (xyz_pt.y()-eval_point.y())*(xyz_pt.y()-eval_point.y()) +
02210                      (xyz_pt.z()-eval_point.z())*(xyz_pt.z()-eval_point.z()) );             
02211         if ( fabs(dsquared) < fabs(dsquared_test) ) {
02212           t_min = t;
02213           dsquared_test = dsquared;
02214         }      
02215         t += 0.2;
02216       }
02217       double dderiva, dderivb;
02218       if ( t_min == -1.00 ) {
02219       //  Check whether the slope changed signs between -1.0 and -0.8.
02220       //  If so, the min must have occurred in this interval -- set
02221       //  t_min to -0.9.
02222         t = -1.0;
02223         edge_ptr->evaluate_single(t,eval_point_ptr);      
02224         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02225         dderiva = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02226                        (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02227                        (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02228         t = -0.8;
02229         edge_ptr->evaluate_single(t,eval_point_ptr);      
02230         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02231         dderivb = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02232                        (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02233                        (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02234         if ( dderiva*dderivb < 0.0 ) t_min = -0.9;
02235    
02236       } else if ( t_min == 1.00 ) {
02237       //  Check the other end of the range, similarly.
02238         t = 1.0;
02239         edge_ptr->evaluate_single(t,eval_point_ptr);      
02240         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02241         dderiva = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02242                         (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02243                         (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02244         t = 0.8;
02245         edge_ptr->evaluate_single(t,eval_point_ptr);      
02246         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02247         dderivb = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02248                         (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02249                         (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02250         if ( dderiva*dderivb < 0.0 ) t_min = 0.9;  
02251       }      
02252       t = t_min;
02253       if ( (t > -1.0) && (t < 1.0) ) {
02254         int mm;
02255         mm = 0;
02256         while ( mm < 10 ) {  //  to avoid possible infinite loop
02257           mm++;
02258           edge_ptr->evaluate_single(t,eval_point_ptr);      
02259           edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02260           dderiv = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02261                          (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02262                          (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02263 
02264           edge_ptr->evaluate_2nd_derivative(t, second_d_ptr);            
02265   
02266           dderiv2 = (eval_point.x()-xyz_pt.x())*second_d.x() + 
02267                      eval_tangent.x()*eval_tangent.x() +
02268                     (eval_point.y()-xyz_pt.y())*second_d.y() + 
02269                      eval_tangent.y()*eval_tangent.y() +
02270                     (eval_point.z()-xyz_pt.z())*second_d.z() + 
02271                      eval_tangent.z()*eval_tangent.z();
02272 
02273            t = t - dderiv/dderiv2; // Newton-Raphson
02274 
02275            if ( t < -1.0 ) {
02276              t = -1.0;
02277              edge_ptr->evaluate_single(t,eval_point_ptr);      
02278              edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02279              break;
02280            }
02281            if ( t > 1.0 ) {
02282              t = 1.0;
02283              edge_ptr->evaluate_single(t,eval_point_ptr);      
02284              edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02285              break;
02286            }
02287            if ( fabs(dderiv) < 1.e-11 ) break;
02288         }
02289       } else {  //  At an endpoint of the paramaterization.
02290           edge_ptr->evaluate_single(t,eval_point_ptr);       
02291           edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02292       }      
02293       xyzOnEdge[index++] = eval_point.x();
02294       xyzOnEdge[index++] = eval_point.y();
02295       xyzOnEdge[index++] = eval_point.z();
02296       if ( tangent != NULL ) {  
02297         eval_tangent.normalize(); 
02298         tangent[tindex++] = eval_tangent.x();
02299         tangent[tindex++] = eval_tangent.y();
02300         tangent[tindex++] = eval_tangent.z();             
02301       }      
02302     }  
02303   }  
02304 
02305   //  delete the temp arrays
02306 
02307   for (ii=0; ii<numEdge; ii++)
02308     delete edge_array[ii];
02309   for (ii=0; ii<numVert; ii++)
02310     delete point_array[ii];
02311   delete [] point_array;
02312   delete [] edge_array;
02313 
02314 }
02315                        
02316 
02317 //===========================================================================
02318 //  Function: projBezierEdge
02319 //  Purpose:  Project a point to an edge.
02320 //            The routine only finds one point.  Rarely there might be
02321 //            two or more closest points.
02322 //  Date:     01/20/2004             
02323 //  Author:   jfowler
02324 //===========================================================================
02325 void projToBezierEdge( int numEdge, int numVert, int* edgeVert, 
02326                               double* vert, double* edgeCtrlPts,
02327                               int numLocs, double* xyz, 
02328                               double* xyzOnEdge, double* tangent, double* t_value )
02329 {
02330 CubitPoint **point_array = new CubitPoint * [numVert];
02331 CubitFacetEdge **edge_array = new CubitFacetEdge * [numEdge];
02332 CubitFacet **facet_array = NULL;
02333 CubitQuadFacet **qfacet_array = NULL;
02334 DLIList<FacetEntity *> facet_list;
02335 int ii, index, numTri, numQuad;
02336 int *triEdge, *quadEdge;
02337   triEdge = quadEdge = NULL;
02338   numTri = numQuad = 0;
02339   CubitStatus stat = CUBIT_SUCCESS;
02340   stat = build_facets(numTri, numQuad, numEdge, numVert, triEdge, quadEdge, 
02341                       edgeVert, vert, point_array, edge_array, facet_array,
02342                       qfacet_array, facet_list);
02343   if (CUBIT_SUCCESS != stat) {
02344     PRINT_ERROR("Failed to build facets.\n");
02345     return;
02346   }
02347 
02348   // set the control points on edges
02349 
02350   CubitFacetEdge *edge_ptr;
02351     //CubitVector T0, T1, P0, P1;
02352     //CubitVector Pi[3];
02353 
02354   for (ii=0; ii<numEdge; ii++)
02355   {
02356     edge_ptr = edge_array[ii];
02357     index = 3 * (ii * NUM_EDGE_CPTS);
02358     edge_ptr->set_control_points( &(edgeCtrlPts[index]) );
02359   } 
02360 
02361   // do the projections
02362   int jj, kk, xyzindex, tindex;
02363   index = tindex = 0;
02364   double dsquared, dsquared_test, dderiv, dderiv2, t, t_min = 0.0;
02365   CubitVector xyz_pt;//, second2142_d;
02366   t = 0.0;
02367   dsquared_test =  CUBIT_DBL_MAX;  
02368   for ( ii = 0; ii < numEdge; ii++ ) {
02369     edge_ptr = edge_array[ii];
02370     for ( jj = 0; jj < numLocs; jj++ ) {
02371       xyzindex = 3*jj;
02372       xyz_pt.x(xyz[xyzindex]);
02373       xyz_pt.y(xyz[xyzindex+1]);
02374       xyz_pt.z(xyz[xyzindex+2]);    
02375       t = -1.0;
02376       dsquared_test =  CUBIT_DBL_MAX;  
02377       //  Get square of distance from xyz to curve at intervals of 0.2
02378       //  in the parameter t.  Save the closest distance found as a 
02379       //  starting point for Newton-Raphson iterations.  It was found 
02380       //  that an interval of 0.5 was too coarse and that the N-R
02381       //  sometimes converged on a false value.
02382       CubitVector eval_point;
02383       CubitVector *eval_point_ptr = &eval_point;
02384       CubitVector eval_tangent;
02385       CubitVector *eval_tangent_ptr = &eval_tangent;
02386       CubitVector second_d;
02387       CubitVector *second_d_ptr = &second_d;
02388       for ( kk = 0; kk < 11; kk++ ) {
02389         edge_ptr->evaluate_single(t,eval_point_ptr);      
02390         dsquared = ( (xyz_pt.x()-eval_point.x())*(xyz_pt.x()-eval_point.x()) +
02391                      (xyz_pt.y()-eval_point.y())*(xyz_pt.y()-eval_point.y()) +
02392                      (xyz_pt.z()-eval_point.z())*(xyz_pt.z()-eval_point.z()) );             
02393         if ( fabs(dsquared) < fabs(dsquared_test) ) {
02394           t_min = t;
02395           dsquared_test = dsquared;
02396         }      
02397         t += 0.2;
02398       }
02399       double dderiva, dderivb;
02400       if ( t_min == -1.00 ) {
02401       //  Check whether the slope changed signs between -1.0 and -0.8.
02402       //  If so, the min must have occurred in this interval -- set
02403       //  t_min to -0.9.
02404         t = -1.0;
02405         edge_ptr->evaluate_single(t,eval_point_ptr);      
02406         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02407         dderiva = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02408                        (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02409                        (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02410         t = -0.8;
02411         edge_ptr->evaluate_single(t,eval_point_ptr);      
02412         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02413         dderivb = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02414                        (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02415                        (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02416         if ( dderiva*dderivb < 0.0 ) t_min = -0.9;
02417    
02418       } else if ( t_min == 1.00 ) {
02419       //  Check the other end of the range, similarly.
02420         t = 1.0;
02421         edge_ptr->evaluate_single(t,eval_point_ptr);      
02422         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02423         dderiva = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02424                         (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02425                         (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02426         t = 0.8;
02427         edge_ptr->evaluate_single(t,eval_point_ptr);      
02428         edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02429         dderivb = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02430                         (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02431                         (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02432         if ( dderiva*dderivb < 0.0 ) t_min = 0.9;  
02433       }      
02434       t = t_min;
02435       if ( (t > -1.0) && (t < 1.0) ) {
02436         int mm;
02437         mm = 0;
02438         while ( mm < 10 ) {  //  to avoid possible infinite loop
02439           mm++;
02440           edge_ptr->evaluate_single(t,eval_point_ptr);      
02441           edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02442           dderiv = 2.0*( (eval_point.x()-xyz_pt.x())*eval_tangent.x() +
02443                          (eval_point.y()-xyz_pt.y())*eval_tangent.y() +
02444                          (eval_point.z()-xyz_pt.z())*eval_tangent.z() );
02445 
02446           edge_ptr->evaluate_2nd_derivative(t, second_d_ptr);            
02447   
02448           dderiv2 = (eval_point.x()-xyz_pt.x())*second_d.x() + 
02449                      eval_tangent.x()*eval_tangent.x() +
02450                     (eval_point.y()-xyz_pt.y())*second_d.y() + 
02451                      eval_tangent.y()*eval_tangent.y() +
02452                     (eval_point.z()-xyz_pt.z())*second_d.z() + 
02453                      eval_tangent.z()*eval_tangent.z();
02454 
02455            t = t - dderiv/dderiv2; // Newton-Raphson
02456 
02457            if ( t < -1.0 ) {
02458              t = -1.0;
02459              edge_ptr->evaluate_single(t,eval_point_ptr);      
02460              edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02461              break;
02462            }
02463            if ( t > 1.0 ) {
02464              t = 1.0;
02465              edge_ptr->evaluate_single(t,eval_point_ptr);      
02466              edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02467              break;
02468            }
02469            if ( fabs(dderiv) < 1.e-11 ) break;
02470         }
02471       } else {  //  At an endpoint of the paramaterization.
02472           edge_ptr->evaluate_single(t,eval_point_ptr);       
02473           edge_ptr->evaluate_single_tangent(t,eval_tangent_ptr);      
02474       }      
02475       xyzOnEdge[index++] = eval_point.x();
02476       xyzOnEdge[index++] = eval_point.y();
02477       xyzOnEdge[index++] = eval_point.z();
02478       if ( tangent != NULL ) {  
02479         eval_tangent.normalize(); 
02480         tangent[tindex++] = eval_tangent.x();
02481         tangent[tindex++] = eval_tangent.y();
02482         tangent[tindex++] = eval_tangent.z();             
02483       }      
02484     }  
02485   } 
02486    
02487   *t_value = t;
02488   
02489   //  delete the temp arrays
02490 
02491   for (ii=0; ii<numEdge; ii++)
02492     delete edge_array[ii];
02493   for (ii=0; ii<numVert; ii++)
02494     delete point_array[ii];
02495   delete [] point_array;
02496   delete [] edge_array;
02497 
02498 }
02499                        
02500 //===========================================================================
02501 //  Function: constructQuadNormalsFromFile
02502 //  Purpose:  Similar to constructQuadNormals, except that instead of
02503 //            approximating the surface normals and tangents from the input
02504 //            facets, the normals and derivatives are extracted from a file
02505 //            that was exported from Cubit along with the exodus file.  The
02506 //            purpose for this is that Cubit will have the exact surface
02507 //            normals and tangents, instead of the approximations from the
02508 //            geometry.
02509 //  Date:     08/31/2004
02510 //  Author:   mlstate
02511 //===========================================================================
02512 CubitStatus constructQuadNormalsFromFile( const char *filename,
02513                                          int numQuad,
02514                                          int numEdge,
02515                                          int* quadEdge,
02516                                          int* edgeVert,
02517                                          double* quadVertNorms,
02518                                          double* edgeVertTang )
02519 {
02520     CubitStatus status = CUBIT_SUCCESS;
02521     int num_file_quads = 0,
02522         num_file_edges = 0,
02523         *edge_nodes_file = NULL,
02524         *quad_nodes_file = NULL;
02525     double *edge_node_tangs_file = NULL,
02526            *quad_node_norms_file = NULL;
02527 
02528     status = readMBGNormalFile( filename,
02529                                 NULL,
02530                                 &num_file_quads,
02531                                 &num_file_edges,
02532                                 &edge_nodes_file,
02533                                 &edge_node_tangs_file,
02534                                 NULL,
02535                                 NULL,
02536                                 &quad_nodes_file,
02537                                 &quad_node_norms_file );
02538 
02539     if ( status == CUBIT_SUCCESS && numEdge > 0 )
02540     {
02541         status = extractEdgeTangsFromFile( num_file_edges,
02542                                            edge_nodes_file,
02543                                            edge_node_tangs_file,
02544                                            numEdge,
02545                                            edgeVert,
02546                                            edgeVertTang );
02547     }
02548     if ( status == CUBIT_SUCCESS && numQuad > 0 )
02549     {
02550         status = extractQuadNormalsFromFile( num_file_quads,
02551                                              quad_nodes_file,
02552                                              quad_node_norms_file,
02553                                              numQuad,
02554                                              edgeVert,
02555                                              quadEdge,
02556                                              quadVertNorms );
02557     }
02558 
02559     delete [] edge_nodes_file;
02560     delete [] quad_nodes_file;
02561     delete [] edge_node_tangs_file;
02562     delete [] quad_node_norms_file;
02563     return status;
02564 }
02565                        
02566 //===========================================================================
02567 //  Function: constructTriNormalsFromFile
02568 //  Purpose:  Similar to constructTriNormals, except that instead of
02569 //            approximating the surface normals and tangents from the input
02570 //            facets, the normals and derivatives are extracted from a file
02571 //            that was exported from Cubit along with the exodus file.  The
02572 //            purpose for this is that Cubit will have the exact surface
02573 //            normals and tangents, instead of the approximations from the
02574 //            geometry.
02575 //  Date:     08/31/2004
02576 //  Author:   mlstate
02577 //===========================================================================
02578 CubitStatus constructTriNormalsFromFile( const char *filename,
02579                                          int numTri,
02580                                          int numEdge,
02581                                          int* triEdge,
02582                                          int* edgeVert,
02583                                          double* triVertNorms,
02584                                          double* edgeVertTang )
02585 {
02586     CubitStatus status = CUBIT_SUCCESS;
02587     int num_file_tris = 0,
02588         num_file_edges = 0,
02589         *edge_nodes_file = NULL,
02590         *tri_nodes_file = NULL;
02591     double *edge_node_tangs_file = NULL,
02592            *tri_node_norms_file = NULL;
02593 
02594     status = readMBGNormalFile( filename,
02595                                 &num_file_tris,
02596                                 NULL,
02597                                 &num_file_edges,
02598                                 &edge_nodes_file,
02599                                 &edge_node_tangs_file,
02600                                 &tri_nodes_file,
02601                                 &tri_node_norms_file,
02602                                 NULL,
02603                                 NULL );
02604 
02605     if ( status == CUBIT_SUCCESS && numEdge > 0 )
02606     {
02607         status = extractEdgeTangsFromFile( num_file_edges,
02608                                            edge_nodes_file,
02609                                            edge_node_tangs_file,
02610                                            numEdge,
02611                                            edgeVert,
02612                                            edgeVertTang );
02613     }
02614     if ( status == CUBIT_SUCCESS && numTri > 0 )
02615     {
02616         status = extractTriNormalsFromFile( num_file_tris,
02617                                             tri_nodes_file,
02618                                             tri_node_norms_file,
02619                                             numTri,
02620                                             edgeVert,
02621                                             triEdge,
02622                                             triVertNorms );
02623     }
02624 
02625     delete [] edge_nodes_file;
02626     delete [] tri_nodes_file;
02627     delete [] edge_node_tangs_file;
02628     delete [] tri_node_norms_file;
02629     return status;
02630 }
02631 
02632 //===========================================================================
02633 //  Function: extractEdgeTangsFromFile
02634 //  Purpose:  extract from the file the edge tangents.
02635 //  Date:     08/31/2004
02636 //  Author:   mlstate
02637 //===========================================================================
02638 static CubitStatus extractEdgeTangsFromFile
02639 (
02640     int num_file_edges,
02641     int *edge_nodes_file,
02642     double *edge_node_tangs_file,
02643     int numEdge,
02644     int *edgeVert,
02645     double *edgeVertTang
02646 )
02647 {
02648 #define MIN_EDGE_NODE_ID( n1, n2 ) n1 < n2 ? n1 : n2
02649     CubitStatus status = CUBIT_SUCCESS;
02650     DLIList<int *> *edge_hash = new DLIList<int*> [num_file_edges];
02651     int * edge_data = new int [num_file_edges * 3 ];
02652 
02653     int i;
02654     
02655     for ( i = 0; i < num_file_edges; i++ )
02656     {
02657         int *data = &(edge_data[i*3]);
02658 
02659         data[0] = edge_nodes_file[i*2];
02660         data[1] = edge_nodes_file[i*2 + 1];
02661         data[2] = i;
02662 
02663         int key = MIN_EDGE_NODE_ID( data[0], data[1] );
02664         key %= num_file_edges;
02665         edge_hash[key].append( data );
02666     }
02667 
02668     for ( i = 0; i < numEdge; i++ )
02669     {
02670         CubitBoolean found = CUBIT_FALSE;
02671         int *this_edge = &(edgeVert[i*2]);
02672         int key = MIN_EDGE_NODE_ID( this_edge[0], this_edge[1] );
02673         key %= num_file_edges;
02674         for ( int j = 0; j < edge_hash[key].size(); j++ )
02675         {
02676             int *data = edge_hash[key].get_and_step();
02677             if ( data[0] == this_edge[0] &&
02678                  data[1] == this_edge[1] )
02679             {
02680                 edgeVertTang[ i*2*3     ] = edge_node_tangs_file[ data[2] * 2 * 3 ];
02681                 edgeVertTang[ i*2*3 + 1 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 1 ];
02682                 edgeVertTang[ i*2*3 + 2 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 2 ];
02683                 edgeVertTang[ i*2*3 + 3 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 3 ];
02684                 edgeVertTang[ i*2*3 + 4 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 4 ];
02685                 edgeVertTang[ i*2*3 + 5 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 5 ];
02686                 found = CUBIT_TRUE;
02687                 break;
02688             }
02689             else if ( data[0] == this_edge[1] &&
02690                       data[1] == this_edge[0] )
02691             {
02692                 edgeVertTang[ i*2*3     ] = edge_node_tangs_file[ data[2] * 2 * 3 + 3 ];
02693                 edgeVertTang[ i*2*3 + 1 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 4 ];
02694                 edgeVertTang[ i*2*3 + 2 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 5 ];
02695                 edgeVertTang[ i*2*3 + 3 ] = edge_node_tangs_file[ data[2] * 2 * 3  ];
02696                 edgeVertTang[ i*2*3 + 4 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 1 ];
02697                 edgeVertTang[ i*2*3 + 5 ] = edge_node_tangs_file[ data[2] * 2 * 3 + 2 ];
02698                 found = CUBIT_TRUE;
02699                 break;
02700             }
02701         }
02702         if ( found == CUBIT_FALSE )
02703         {
02704             status = CUBIT_FAILURE;
02705             break;
02706         }
02707     }
02708     delete [] edge_data;
02709     delete [] edge_hash;
02710 
02711     return status;
02712  }
02713 
02714 //===========================================================================
02715 //  Function: extractTriNormalsFromFile
02716 //  Purpose:  extract from the file the normals for a block of tri elements.
02717 //  Date:     08/31/2004
02718 //  Author:   mlstate
02719 //===========================================================================
02720 static CubitStatus extractTriNormalsFromFile
02721 (
02722     int num_file_tris,
02723     int *tri_nodes_file,
02724     double *tri_node_norms_file,
02725     int numTri,
02726     int *edgeVert,
02727     int *triEdge,
02728     double *triVertNorms
02729 )
02730 {
02731 #define MIN_TRI_NODE_ID( n1, n2, n3 ) n1 < n2 ?  n1 < n3 ? n1 : n3 : n2 < n3 ? n2 : n3
02732 
02733     CubitStatus status = CUBIT_SUCCESS;
02734     DLIList< int* > *tri_hash = new DLIList< int*>[num_file_tris];
02735     int * tri_data = new int [num_file_tris * 4 ];
02736 
02737     int i;
02738     
02739     for ( i = 0; i < num_file_tris; i++ )
02740     {
02741         int *data = &(tri_data[i*4]);
02742 
02743         data[0] = tri_nodes_file[i*3];
02744         data[1] = tri_nodes_file[i*3 + 1];
02745         data[2] = tri_nodes_file[i*3 + 2];
02746         data[3] = i;
02747 
02748         int key = MIN_TRI_NODE_ID( data[0], data[1], data[2] );
02749         key %= num_file_tris;
02750         tri_hash[key].append( data );
02751     }
02752 
02753     for ( i = 0; i < numTri; i++ )
02754     {
02755         CubitBoolean found = CUBIT_FALSE;
02756         int this_tri[3];
02757         
02758         constructTriVerts( &triEdge[i*3], edgeVert, this_tri );
02759 
02760         int key = MIN_TRI_NODE_ID( this_tri[0], this_tri[1], this_tri[2] );
02761         key %= num_file_tris;
02762         for ( int j = 0; j < tri_hash[key].size(); j++ )
02763         {
02764             int *data = tri_hash[key].get_and_step();
02765 
02766             if ( this_tri[0] == data[0] &&
02767                  this_tri[1] == data[1] &&
02768                  this_tri[2] == data[2] )
02769             {
02770                 triVertNorms[ i*3*3     ] = tri_node_norms_file[ data[3]*3*3     ];
02771                 triVertNorms[ i*3*3 + 1 ] = tri_node_norms_file[ data[3]*3*3 + 1 ];
02772                 triVertNorms[ i*3*3 + 2 ] = tri_node_norms_file[ data[3]*3*3 + 2 ];
02773                 triVertNorms[ i*3*3 + 3 ] = tri_node_norms_file[ data[3]*3*3 + 3 ];
02774                 triVertNorms[ i*3*3 + 4 ] = tri_node_norms_file[ data[3]*3*3 + 4 ];
02775                 triVertNorms[ i*3*3 + 5 ] = tri_node_norms_file[ data[3]*3*3 + 5 ];
02776                 triVertNorms[ i*3*3 + 6 ] = tri_node_norms_file[ data[3]*3*3 + 6 ];
02777                 triVertNorms[ i*3*3 + 7 ] = tri_node_norms_file[ data[3]*3*3 + 7 ];
02778                 triVertNorms[ i*3*3 + 8 ] = tri_node_norms_file[ data[3]*3*3 + 8 ];
02779                 found = CUBIT_TRUE;
02780                 break;
02781             }
02782             else if ( this_tri[0] == data[1] &&
02783                       this_tri[1] == data[2] &&
02784                       this_tri[2] == data[0] )
02785             {
02786                 triVertNorms[ i*3*3     ] = tri_node_norms_file[ data[3]*3*3 + 3 ];
02787                 triVertNorms[ i*3*3 + 1 ] = tri_node_norms_file[ data[3]*3*3 + 4 ];
02788                 triVertNorms[ i*3*3 + 2 ] = tri_node_norms_file[ data[3]*3*3 + 5 ];
02789                 triVertNorms[ i*3*3 + 3 ] = tri_node_norms_file[ data[3]*3*3 + 6 ];
02790                 triVertNorms[ i*3*3 + 4 ] = tri_node_norms_file[ data[3]*3*3 + 7 ];
02791                 triVertNorms[ i*3*3 + 5 ] = tri_node_norms_file[ data[3]*3*3 + 8 ];
02792                 triVertNorms[ i*3*3 + 6 ] = tri_node_norms_file[ data[3]*3*3     ];
02793                 triVertNorms[ i*3*3 + 7 ] = tri_node_norms_file[ data[3]*3*3 + 1 ];
02794                 triVertNorms[ i*3*3 + 8 ] = tri_node_norms_file[ data[3]*3*3 + 2 ];
02795                 found = CUBIT_TRUE;
02796                 break;
02797             }
02798             else if ( this_tri[0] == data[2] &&
02799                       this_tri[1] == data[0] &&
02800                       this_tri[2] == data[1] )
02801             {
02802                 triVertNorms[ i*3*3     ] = tri_node_norms_file[ data[3]*3*3 + 6 ];
02803                 triVertNorms[ i*3*3 + 1 ] = tri_node_norms_file[ data[3]*3*3 + 7 ];
02804                 triVertNorms[ i*3*3 + 2 ] = tri_node_norms_file[ data[3]*3*3 + 8 ];
02805                 triVertNorms[ i*3*3 + 3 ] = tri_node_norms_file[ data[3]*3*3     ];
02806                 triVertNorms[ i*3*3 + 4 ] = tri_node_norms_file[ data[3]*3*3 + 1 ];
02807                 triVertNorms[ i*3*3 + 5 ] = tri_node_norms_file[ data[3]*3*3 + 2 ];
02808                 triVertNorms[ i*3*3 + 6 ] = tri_node_norms_file[ data[3]*3*3 + 3 ];
02809                 triVertNorms[ i*3*3 + 7 ] = tri_node_norms_file[ data[3]*3*3 + 4 ];
02810                 triVertNorms[ i*3*3 + 8 ] = tri_node_norms_file[ data[3]*3*3 + 5 ];
02811                 found = CUBIT_TRUE;
02812                 break;
02813             }
02814         }
02815         if ( found == CUBIT_FALSE )
02816         {
02817             status = CUBIT_FAILURE;
02818             break;
02819         }
02820     }
02821 
02822     delete [] tri_data;
02823     delete [] tri_hash;
02824     return status;
02825 }
02826 
02827 //===========================================================================
02828 //  Function: extractQuadNormalsFromFile
02829 //  Purpose:  extract from the file the normals for a block of quad elements.
02830 //  Date:     08/31/2004
02831 //  Author:   mlstate
02832 //===========================================================================
02833 static CubitStatus extractQuadNormalsFromFile
02834 (
02835     int num_file_quads,
02836     int *quad_nodes_file,
02837     double *quad_node_norms_file,
02838     int numQuad,
02839     int *edgeVert,
02840     int *quadEdge,
02841     double *quadVertNorms
02842 )
02843 {
02844 #define MIN_QUAD_NODE_ID( n1, n2, n3, n4 ) MIN_EDGE_NODE_ID( (MIN_EDGE_NODE_ID( n1, n2 )), (MIN_EDGE_NODE_ID( n3, n4 )) )
02845 
02846     CubitStatus status = CUBIT_SUCCESS;
02847     DLIList< int* > *quad_hash = new DLIList< int*>[num_file_quads];
02848     int * quad_data = new int [num_file_quads * 5 ];
02849 
02850     int i;
02851     for ( i = 0; i < num_file_quads; i++ )
02852     {
02853         int *data = &(quad_data[i*5]);
02854 
02855         data[0] = quad_nodes_file[i*4];
02856         data[1] = quad_nodes_file[i*4 + 1];
02857         data[2] = quad_nodes_file[i*4 + 2];
02858         data[3] = quad_nodes_file[i*4 + 3];
02859         data[4] = i;
02860 
02861         int key = MIN_QUAD_NODE_ID( data[0], data[1], data[2], data[3] );
02862         key %= num_file_quads;
02863         quad_hash[key].append( data );
02864     }
02865 
02866     for ( i = 0; i < numQuad; i++ )
02867     {
02868         CubitBoolean found = CUBIT_FALSE;
02869         int this_quad[4];
02870         
02871         constructQuadVerts( &quadEdge[i*4], edgeVert, this_quad );
02872 
02873         int key = MIN_QUAD_NODE_ID( this_quad[0], this_quad[1], this_quad[2], this_quad[3] );
02874         key %= num_file_quads;
02875         for ( int j = 0; j < quad_hash[key].size(); j++ )
02876         {
02877             int *data = quad_hash[key].get_and_step();
02878 
02879             if ( this_quad[0] == data[0] &&
02880                  this_quad[1] == data[1] &&
02881                  this_quad[2] == data[2] &&
02882                  this_quad[3] == data[3] )
02883             {
02884                 quadVertNorms[ i*4*3     ] = quad_node_norms_file[ data[4]*4*3     ];
02885                 quadVertNorms[ i*4*3 + 1 ] = quad_node_norms_file[ data[4]*4*3 + 1 ];
02886                 quadVertNorms[ i*4*3 + 2 ] = quad_node_norms_file[ data[4]*4*3 + 2 ];
02887                 quadVertNorms[ i*4*3 + 3 ] = quad_node_norms_file[ data[4]*4*3 + 3 ];
02888                 quadVertNorms[ i*4*3 + 4 ] = quad_node_norms_file[ data[4]*4*3 + 4 ];
02889                 quadVertNorms[ i*4*3 + 5 ] = quad_node_norms_file[ data[4]*4*3 + 5 ];
02890                 quadVertNorms[ i*4*3 + 6 ] = quad_node_norms_file[ data[4]*4*3 + 6 ];
02891                 quadVertNorms[ i*4*3 + 7 ] = quad_node_norms_file[ data[4]*4*3 + 7 ];
02892                 quadVertNorms[ i*4*3 + 8 ] = quad_node_norms_file[ data[4]*4*3 + 8 ];
02893                 quadVertNorms[ i*4*3 + 9 ] = quad_node_norms_file[ data[4]*4*3 + 9 ];
02894                 quadVertNorms[ i*4*3 + 10] = quad_node_norms_file[ data[4]*4*3 + 10];
02895                 quadVertNorms[ i*4*3 + 11] = quad_node_norms_file[ data[4]*4*3 + 11];
02896                 found = CUBIT_TRUE;
02897                 break;
02898             }
02899             else if ( this_quad[0] == data[1] &&
02900                       this_quad[1] == data[2] &&
02901                       this_quad[2] == data[3] &&
02902                       this_quad[3] == data[0] )
02903             {
02904                 quadVertNorms[ i*4*3     ] = quad_node_norms_file[ data[4]*4*3 + 3 ];
02905                 quadVertNorms[ i*4*3 + 1 ] = quad_node_norms_file[ data[4]*4*3 + 4 ];
02906                 quadVertNorms[ i*4*3 + 2 ] = quad_node_norms_file[ data[4]*4*3 + 5 ];
02907                 quadVertNorms[ i*4*3 + 3 ] = quad_node_norms_file[ data[4]*4*3 + 6 ];
02908                 quadVertNorms[ i*4*3 + 4 ] = quad_node_norms_file[ data[4]*4*3 + 7 ];
02909                 quadVertNorms[ i*4*3 + 5 ] = quad_node_norms_file[ data[4]*4*3 + 8 ];
02910                 quadVertNorms[ i*4*3 + 6 ] = quad_node_norms_file[ data[4]*4*3 + 9 ];
02911                 quadVertNorms[ i*4*3 + 7 ] = quad_node_norms_file[ data[4]*4*3 + 10 ];
02912                 quadVertNorms[ i*4*3 + 8 ] = quad_node_norms_file[ data[4]*4*3 + 11 ];
02913                 quadVertNorms[ i*4*3 + 9 ] = quad_node_norms_file[ data[4]*4*3 + 0 ];
02914                 quadVertNorms[ i*4*3 + 10] = quad_node_norms_file[ data[4]*4*3 + 1 ];
02915                 quadVertNorms[ i*4*3 + 11] = quad_node_norms_file[ data[4]*4*3 + 2 ];
02916                 found = CUBIT_TRUE;
02917                 break;
02918             }
02919             else if ( this_quad[0] == data[2] &&
02920                       this_quad[1] == data[3] &&
02921                       this_quad[2] == data[0] &&
02922                       this_quad[3] == data[1] )
02923             {
02924                 quadVertNorms[ i*4*3     ] = quad_node_norms_file[ data[4]*4*3 + 6 ];
02925                 quadVertNorms[ i*4*3 + 1 ] = quad_node_norms_file[ data[4]*4*3 + 7 ];
02926                 quadVertNorms[ i*4*3 + 2 ] = quad_node_norms_file[ data[4]*4*3 + 8 ];
02927                 quadVertNorms[ i*4*3 + 3 ] = quad_node_norms_file[ data[4]*4*3 + 9 ];
02928                 quadVertNorms[ i*4*3 + 4 ] = quad_node_norms_file[ data[4]*4*3 + 10 ];
02929                 quadVertNorms[ i*4*3 + 5 ] = quad_node_norms_file[ data[4]*4*3 + 11 ];
02930                 quadVertNorms[ i*4*3 + 6 ] = quad_node_norms_file[ data[4]*4*3 + 0 ];
02931                 quadVertNorms[ i*4*3 + 7 ] = quad_node_norms_file[ data[4]*4*3 + 1 ];
02932                 quadVertNorms[ i*4*3 + 8 ] = quad_node_norms_file[ data[4]*4*3 + 2 ];
02933                 quadVertNorms[ i*4*3 + 9 ] = quad_node_norms_file[ data[4]*4*3 + 3 ];
02934                 quadVertNorms[ i*4*3 + 10] = quad_node_norms_file[ data[4]*4*3 + 4 ];
02935                 quadVertNorms[ i*4*3 + 11] = quad_node_norms_file[ data[4]*4*3 + 5 ];
02936                 found = CUBIT_TRUE;
02937                 break;
02938             }
02939             else if (  this_quad[0] == data[3] &&
02940                        this_quad[1] == data[0] &&
02941                        this_quad[2] == data[1] &&
02942                        this_quad[3] == data[2] )
02943             {
02944                 quadVertNorms[ i*4*3     ] = quad_node_norms_file[ data[4]*4*3 + 9 ];
02945                 quadVertNorms[ i*4*3 + 1 ] = quad_node_norms_file[ data[4]*4*3 + 10 ];
02946                 quadVertNorms[ i*4*3 + 2 ] = quad_node_norms_file[ data[4]*4*3 + 11 ];
02947                 quadVertNorms[ i*4*3 + 3 ] = quad_node_norms_file[ data[4]*4*3 + 0 ];
02948                 quadVertNorms[ i*4*3 + 4 ] = quad_node_norms_file[ data[4]*4*3 + 1 ];
02949                 quadVertNorms[ i*4*3 + 5 ] = quad_node_norms_file[ data[4]*4*3 + 2 ];
02950                 quadVertNorms[ i*4*3 + 6 ] = quad_node_norms_file[ data[4]*4*3 + 3 ];
02951                 quadVertNorms[ i*4*3 + 7 ] = quad_node_norms_file[ data[4]*4*3 + 4 ];
02952                 quadVertNorms[ i*4*3 + 8 ] = quad_node_norms_file[ data[4]*4*3 + 5 ];
02953                 quadVertNorms[ i*4*3 + 9 ] = quad_node_norms_file[ data[4]*4*3 + 6 ];
02954                 quadVertNorms[ i*4*3 + 10] = quad_node_norms_file[ data[4]*4*3 + 7 ];
02955                 quadVertNorms[ i*4*3 + 11] = quad_node_norms_file[ data[4]*4*3 + 8 ];
02956                 found = CUBIT_TRUE;
02957                 break;
02958             }
02959         }
02960         if ( found == CUBIT_FALSE )
02961         {
02962             status = CUBIT_FAILURE;
02963             break;
02964         }
02965     }
02966 
02967     delete [] quad_data;
02968     delete [] quad_hash;
02969     return status;
02970 }
02971 
02972 //===========================================================================
02973 //  Function: constructTriVerts
02974 //  Purpose:  Given an array of edge indices for a tri element and an array
02975 //            of edge to vertex info, construct an array containing the vertices
02976 //            on the input tri.
02977 //  Date:     08/31/2004
02978 //  Author:   mlstate
02979 //===========================================================================
02980 static void constructTriVerts
02981 (
02982     int triEdge[3],  // <I> Array of edges on this tri.
02983     int *edgeVert,   // <I> Array of edge to vertex info.
02984     int triVert[3]   // <O> The vertices on the input tri.
02985 )
02986 {
02987     int *verts_edge1 = &(edgeVert[triEdge[0]*2]);
02988     int *verts_edge2 = &(edgeVert[triEdge[1]*2]);
02989     int *verts_edge3 = &(edgeVert[triEdge[2]*2]);
02990 
02991     if ( verts_edge1[0] != verts_edge2[0] &&
02992          verts_edge1[0] != verts_edge2[1] )
02993     {
02994         triVert[0] = verts_edge1[0];
02995         triVert[1] = verts_edge1[1];
02996     }
02997     else
02998     {
02999         triVert[0] = verts_edge1[1];
03000         triVert[1] = verts_edge1[0];
03001     }
03002     if ( verts_edge3[0] != verts_edge2[0] &&
03003          verts_edge3[0] != verts_edge2[1] )
03004     {
03005         triVert[2] = verts_edge3[1];
03006     }
03007     else
03008     {
03009         triVert[2] = verts_edge3[0];
03010     }
03011 }
03012 
03013 //===========================================================================
03014 //  Function: constructQuadVerts
03015 //  Purpose:  Given an array of edge indices for a quad element and an array
03016 //            of edge to vertex info, construct an array containing the vertices
03017 //            on the input quad.
03018 //  Date:     08/31/2004
03019 //  Author:   mlstate
03020 //===========================================================================
03021 static void constructQuadVerts
03022 (
03023     int QuadEdge[4],  // <I> Array of edges on this quad.
03024     int *edgeVert,    // <I> Array of edge to vertex info.
03025     int QuadVert[4]   // <O> The vertices on the input quad.
03026 )
03027 {
03028     int *verts_edge1 = &(edgeVert[QuadEdge[0]*2]);
03029     int *verts_edge2 = &(edgeVert[QuadEdge[1]*2]);
03030     int *verts_edge3 = &(edgeVert[QuadEdge[2]*2]);
03031     int *verts_edge4 = &(edgeVert[QuadEdge[3]*2]);
03032 
03033     if ( verts_edge1[0] != verts_edge2[0] &&
03034          verts_edge1[0] != verts_edge2[1] )
03035     {
03036         QuadVert[0] = verts_edge1[0];
03037         QuadVert[1] = verts_edge1[1];
03038     }
03039     else
03040     {
03041         QuadVert[0] = verts_edge1[1];
03042         QuadVert[1] = verts_edge1[0];
03043     }
03044 
03045     if ( verts_edge3[0] != verts_edge4[0] &&
03046          verts_edge3[0] != verts_edge4[1] )
03047     {
03048         QuadVert[2] = verts_edge3[0];
03049         QuadVert[3] = verts_edge3[1];
03050     }
03051     else
03052     {
03053         QuadVert[2] = verts_edge3[1];
03054         QuadVert[3] = verts_edge3[0];
03055     }
03056 }
03057 
03058 //===========================================================================
03059 //  Function: readMBGNormalFile
03060 //  Purpose:  Read the raw data out of the Normal/Tangents file.
03061 //  Date:     08/31/2004
03062 //  Author:   mlstate
03063 //===========================================================================
03064 CubitStatus readMBGNormalFile
03065 (
03066     const char *filename,
03067     int *num_tris,                // <O>
03068     int *num_quads,               // <O>
03069     int *num_edges,               // <O>
03070     int **edge_nodes,             // <OF> len = 2 * num_edges
03071     double **edge_node_tangents,  // <OF> len = 3 * 2 * num_edges
03072     int **tri_nodes,              // <OF> len = 3 * num_tris
03073     double **tri_node_normals,    // <OF> len = 3 * 3 * num_tris
03074     int **quad_nodes,             // <OF> len = 4 * num_quads
03075     double **quad_node_normals    // <OF> len = 3 * 4 * num_quads
03076 )
03077 {
03078 #define MAX_FILE_LINE 512
03079 
03080     int num_t = 0;
03081     int num_q = 0;
03082     int curr_quad = 0;
03083     int curr_tri = 0;
03084     int num_e = 0;
03085     int n[4];
03086 
03087     if ( num_tris           ) *num_tris = 0;
03088     if ( num_quads          ) *num_quads = 0;
03089     if ( num_edges          ) *num_edges = 0;
03090     if ( edge_nodes         ) *edge_nodes = NULL;
03091     if ( edge_node_tangents ) *edge_node_tangents = NULL;
03092     if ( tri_nodes          ) *tri_nodes = NULL;
03093     if ( quad_nodes         ) *quad_nodes = NULL;
03094     if ( tri_node_normals   ) *tri_node_normals = NULL;
03095     if ( quad_node_normals  ) *quad_node_normals = NULL;
03096 
03097     if ( strlen( filename ) == 0 )
03098     {
03099         PRINT_ERROR(" No filename specified\n");
03100         return CUBIT_FAILURE;
03101     }
03102 
03103     ifstream geom_file( filename );
03104 
03105     char fileline[MAX_FILE_LINE];
03106 
03107     while ( geom_file.getline( fileline, MAX_FILE_LINE ) )
03108     {
03109         int block_id,
03110             num_corners,
03111             num_elems;
03112 
03113         if ( !strcmp( fileline, END_OF_BLOCKS ) ) break;
03114         if ( 3 != sscanf( fileline, "BLK %d %d %d",
03115                           &block_id, &num_corners, &num_elems ) )
03116         {
03117             PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03118             geom_file.close();
03119             return CUBIT_FAILURE;
03120         }
03121 
03122         if ( num_corners == 3 )
03123         {
03124             checkMemoryAllocations( num_corners, num_elems, &num_t, tri_nodes,
03125                                     tri_node_normals );
03126         }
03127         else if ( num_corners == 4 )
03128         {
03129             checkMemoryAllocations( num_corners, num_elems, &num_q, quad_nodes,
03130                                     quad_node_normals );
03131         }
03132         else
03133         {
03134             PRINT_ERROR( "MBG Normal file can only support quads and tris.\n" );
03135             geom_file.close();
03136             return CUBIT_FAILURE;
03137         }
03138 
03139         for ( int ielem = 0; ielem < num_elems; ielem++ )
03140         {
03141             int *this_nodes = NULL;
03142             double *this_norms = NULL;
03143 
03144             if ( !geom_file.getline( fileline, MAX_FILE_LINE ) )
03145             {
03146                 PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03147                 geom_file.close();
03148                 return CUBIT_FAILURE;
03149             }
03150 
03151             if ( num_corners == 4 )
03152             {
03153                 if ( quad_nodes        ) this_nodes = &((*quad_nodes       )[4     * curr_quad]);
03154                 if ( quad_node_normals ) this_norms = &((*quad_node_normals)[4 * 3 * curr_quad]);
03155                 curr_quad++;
03156             }
03157             else if ( num_corners == 3 )
03158             {
03159                 if ( tri_nodes        ) this_nodes = &((*tri_nodes       )[3     * curr_tri]);
03160                 if ( tri_node_normals ) this_norms = &((*tri_node_normals)[3 * 3 * curr_tri]);
03161                 curr_tri++;
03162             }
03163             if ( ( num_corners == 4 && 4 != sscanf( fileline, "%d %d %d %d",
03164                                                     &n[0], &n[1], &n[2], &n[3]  ) ) ||
03165                  ( num_corners == 3 && 3 != sscanf( fileline, "%d %d %d",
03166                                                     &n[0], &n[1], &n[2] ) ) )
03167             {
03168                 PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03169                 geom_file.close();
03170                 return CUBIT_FAILURE;
03171             }
03172 
03173             if ( this_nodes )
03174                 memcpy( this_nodes, n, num_corners * sizeof( int ) );
03175 
03176             // Read in the surface normals at the vertices of this element.
03177 
03178             for ( int icorner = 0; icorner < num_corners; icorner++ )
03179             {
03180                 if ( !geom_file.getline( fileline, MAX_FILE_LINE ) )
03181                 {
03182                     PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03183                     geom_file.close();
03184                     return CUBIT_FAILURE;
03185                 }
03186                 if ( this_norms )
03187                 {
03188                     double norm[3];
03189 
03190                     if ( 3 != sscanf( fileline, "%lf %lf %lf", &norm[0], &norm[1], &norm[2] ) )
03191                     {
03192                         PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03193                         geom_file.close();
03194                         return CUBIT_FAILURE;
03195                     }
03196                     memcpy( &(this_norms[icorner * 3]), norm, 3 * sizeof( double ) );
03197                 }
03198             }
03199         }
03200     }
03201 
03202     if ( !geom_file.getline( fileline, MAX_FILE_LINE ) )
03203     {
03204         PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03205         geom_file.close();
03206         return CUBIT_FAILURE;
03207     }
03208 
03209     if ( 1 != sscanf( fileline, "%d", &num_e ) )
03210     {
03211         PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03212         geom_file.close();
03213         return CUBIT_FAILURE;
03214     }
03215 
03216     if ( edge_nodes || edge_node_tangents )
03217     {
03218         if ( edge_nodes         ) *edge_nodes         = new int [ 2 * num_e ];
03219         if ( edge_node_tangents ) *edge_node_tangents = new double [3 * 2 * num_e];
03220 
03221         for ( int iedge = 0; iedge < num_e; iedge++ )
03222         {
03223             if ( !geom_file.getline( fileline, MAX_FILE_LINE ) )
03224             {
03225                 PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03226                 geom_file.close();
03227                 return CUBIT_FAILURE;
03228             }
03229             if ( 2 != sscanf( fileline, "%d %d", &n[0], &n[1] ) )
03230             {
03231                 PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03232                 geom_file.close();
03233                 return CUBIT_FAILURE;
03234             }
03235 
03236             for ( int iend = 0; iend < 2; iend++ )
03237             {
03238                 double d[3];
03239 
03240                 if ( edge_nodes )
03241                 {
03242                     (*edge_nodes)[iedge * 2 + iend ] = n[iend];
03243                 }
03244 
03245                 if ( !geom_file.getline( fileline, MAX_FILE_LINE ) )
03246                 {
03247                     PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03248                     geom_file.close();
03249                     return CUBIT_FAILURE;
03250                 }
03251                 if ( 3 != sscanf( fileline, "%lf %lf %lf", &d[0], &d[1], &d[2] ) )
03252                 {
03253                     PRINT_ERROR( "MBG Normal file has incorrect file format\n" );
03254                     geom_file.close();
03255                     return CUBIT_FAILURE;
03256                 }
03257 
03258                 int index = ((iedge*2) + iend) * 3;
03259                 if ( edge_node_tangents )
03260                 {
03261                     (*edge_node_tangents)[ index     ] = d[0];
03262                     (*edge_node_tangents)[ index + 1 ] = d[1];
03263                     (*edge_node_tangents)[ index + 2 ] = d[2];
03264                 }
03265             }
03266         }
03267     }
03268 
03269     if ( num_quads ) *num_quads = num_q;
03270     if ( num_tris  ) *num_tris  = num_t;
03271     if ( num_edges ) *num_edges = num_e;
03272 
03273     geom_file.close();
03274     return CUBIT_SUCCESS;
03275 }
03276 
03277 //===========================================================================
03278 //  Function: checkMemoryAllocations
03279 //  Purpose:  While reading data out of the Normal/Tangents file, reallocate
03280 //            Arrays to ensure that we have enough memory.
03281 //  Date:     08/31/2004
03282 //  Author:   mlstate
03283 //===========================================================================
03284 static void checkMemoryAllocations
03285 (
03286     int num_nodes_per_elem,
03287     int additional_num_elems,
03288     int *num_elems,
03289     int **nodes,
03290     double **node_normals
03291 )
03292 {
03293     int old_num = *num_elems;
03294 
03295     *num_elems += additional_num_elems;
03296     if ( node_normals && *node_normals )
03297     {
03298         double *new_mem = new double [ *num_elems * 3 * num_nodes_per_elem ];
03299         memcpy( new_mem, *node_normals,  old_num * 3 * num_nodes_per_elem * sizeof( double ) );
03300         delete *node_normals;
03301         *node_normals = new_mem;
03302     }
03303     else if ( node_normals )
03304     {
03305         *node_normals = new double [ *num_elems * num_nodes_per_elem * 3 ];
03306     }
03307 
03308     if ( nodes && *nodes )
03309     {
03310         int *new_mem = new int [ *num_elems * num_nodes_per_elem ];
03311         memcpy( new_mem, *nodes,  old_num * num_nodes_per_elem * sizeof( int ) );
03312         delete *nodes;
03313         *nodes = new_mem;
03314     }
03315     else if ( nodes )
03316     {
03317         *nodes = new int [ *num_elems * num_nodes_per_elem ];
03318     }
03319 }
03320 
03321 // This code classify the type of edge with respect to convexity,
03322 // C0 or C1 continuity, and topology.  The returned int uses the following
03323 // number code:
03324 // -1 = non-convex, manifold edge
03325 //  0 = non-feature, manifold  edge (ie, C1 edge)
03326 //  1 = convex, manifold edge
03327 //  2 = non-manifold edge
03328 //  Numbers higher than two are error codes
03329 //  3 = edge only attached to a single facet (boundary of 2-d sheet).
03330 static int classify_local_convexity_at_edge(CubitFacetEdge *edge_ptr)
03331 {
03332     //easy check... if not a feature return 0
03333   if(!edge_ptr->is_feature())
03334     return 0;
03335   
03336     //else we know that we are at feature... determine the type
03337   CubitFacet *adj_facet_1 = NULL, *adj_facet_2 = NULL;
03338   CubitPoint *point_1 = NULL, *point_2 = NULL, *point_3 = NULL;
03339   CubitVector v2, v3;
03340   CubitVector v0 = edge_ptr->point(0)->coordinates();
03341   CubitVector v1 = edge_ptr->point(1)->coordinates();
03342   DLIList<CubitFacet*> adj_facets;
03343   edge_ptr->facets(adj_facets);
03344   
03345     //if non-manifold edge, mark as two
03346   if(adj_facets.size() > 2){
03347     return 2;
03348   }
03349     //if boundary on 2-d (shouldn't happen) mark as 3
03350     //anything above 2 is essentially an error code...
03351   else if(adj_facets.size() < 2){
03352     return 3;
03353   }
03354     //otherwise, determine whether it is 1 or -1
03355     //this else should be a separate function:  classify_edge_convexity()
03356   adj_facet_1 = adj_facets.get_and_step();
03357   int edge_index = adj_facet_1->edge_index(edge_ptr);
03358     //order such that *_1 is the facet that uses
03359     // this edge as the forward edge...
03360   if(adj_facet_1->edge_use(edge_index) == 1)
03361     adj_facet_2 = adj_facets.get();
03362   else{
03363     adj_facet_1 = adj_facets.get();
03364     adj_facets.reset();
03365     adj_facet_2 = adj_facets.get();
03366   }
03367     // order the points so that the Jacobian of the matrix
03368     // will tell use whether we are convex or not at this edge.
03369   adj_facet_1->points(point_1, point_2, point_3);
03370   if(point_1 == edge_ptr->point(0)){
03371     v2 = point_3->coordinates();
03372   }
03373   else if(point_1 == edge_ptr->point(1)){
03374     v2 = point_2->coordinates();
03375   }
03376   else{
03377     v2 = point_1->coordinates();
03378   }
03379   adj_facet_2->points(point_1, point_2, point_3);
03380   if(point_1 == edge_ptr->point(0)){
03381     v3 = point_2->coordinates();
03382   }
03383   else if(point_1 == edge_ptr->point(1)){
03384     v3 = point_3->coordinates();
03385   }
03386   else{
03387     v3 = point_1->coordinates();
03388   }
03389 
03390   v1-=v0;
03391   v2-=v0;
03392   v3-=v0;
03393     //if Jacobian is negative, we are convex here.
03394   if( ( v1 % (v2* v3) ) < 0.0){
03395     return 1;
03396   }
03397   else{
03398     return -1;
03399   }
03400 
03401 }
03402 /*
03403 int test_edges(int num_edge, int* edge_vert, double* my_vert, CubitFacetEdge** edge_array){
03404     int ii;
03405     int return_val = 0.0;
03406   CubitFacetEdgeData *cfed_ptr = NULL;
03407   int my_v1, my_v2;
03408   CubitVector vert1, vert2;
03409   CubitPoint *poi1, *poi2;
03410   CubitVector coord1, coord2;
03411   for(ii=0; ii<num_edge; ++ii){
03412     cfed_ptr = CAST_TO(edge_array[ii], CubitFacetEdgeData);
03413     poi1 = cfed_ptr->start_node();
03414     poi2 = cfed_ptr->end_node();
03415     coord1.set(poi1->x(),poi1->y(),poi1->z());
03416     coord2.set(poi2->x(),poi2->y(),poi2->z());
03417     my_v1 = edge_vert[ii*2];
03418     my_v2 = edge_vert[ii*2+1];
03419     vert1.set(my_vert[my_v1*3],my_vert[my_v1*3+1],my_vert[my_v1*3+2]);
03420     vert2.set(my_vert[my_v2*3],my_vert[my_v2*3+1],my_vert[my_v2*3+2]);
03421     double my_dist = vert1.distance_between(coord1);
03422     if(my_dist > .01){
03423         PRINT_INFO("Edge (%i)(1), distance (%f)\n",ii,my_dist);
03424         return_val++;
03425     }
03426     my_dist = vert2.distance_between(coord2);
03427     if(my_dist > .01){
03428         PRINT_INFO("Edge (%i)(2), distance (%f)\n",ii,my_dist);
03429         return_val++;
03430     }
03431   }
03432   return return_val;
03433 } 
03434 
03435   
03436 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines