cgma
CubitOctree.cpp
Go to the documentation of this file.
00001 #include "CubitOctree.hpp"
00002 #include "RefFace.hpp" 
00003 #include "CubitOctreeConstants.hpp"
00004 
00005 #include "CubitOctreeCell.hpp"
00006 #include "CubitOctreeNode.hpp"  
00007 #include "CubitFacet.hpp" 
00008 #include "CubitFacetEdge.hpp" 
00009 #include "CubitPoint.hpp"
00010 #include "CubitPointData.hpp"
00011 #include "TDOctreeRefFace.hpp"
00012 #include "FacetDataUtil.hpp"
00013 #include "CubitOctreeGenerator.hpp"
00014 #include "OctreeFacetPointData.hpp"
00015 
00016 
00017 
00018 /* -------------------- Methods of CubitOctree ---------------- */
00019 //      Takes the following parameters
00020 //      min_depth : Minimum depth of the CubitOctree (avoids large cubes).
00021 //      max_depht : Maximum depth of the CubitOctree (restricted by memory). 
00022     
00023 CubitOctree::CubitOctree( CubitOctreeGenerator *lattice_gen ){
00024 
00025   CubitOctreeNode::reset_counter();
00026   
00027   cubitOctreeGenerator = lattice_gen;
00028 
00029 
00030   maxSizeGridNode = -CUBIT_DBL_MAX; 
00031   minSizeGridNode = CUBIT_DBL_MAX;
00032 
00033   CubitVector min_corner; 
00034   CubitVector max_corner;  
00035   
00036   lattice_gen->get_bounding_box( min_corner, max_corner );
00037 
00038     // note (vvyas 7/2005): perturbing bbox unsymmetrically to help avoid degenerate cases
00039     // in SAT intersection code
00040   CubitVector min_perturbation(-.0001,-.0002,-.0003);
00041   CubitVector max_perturbation(.0005,.00012,.00028);
00042 
00043   min_corner += min_perturbation;
00044   max_corner += max_perturbation;
00045   
00046   CubitBox bbox(min_corner, max_corner );
00047   
00048   double x_range = bbox.x_range();
00049   double y_range = bbox.y_range();
00050   double z_range = bbox.z_range();
00051   
00052   double max_range;
00053   max_range = ( x_range > y_range ) ? x_range : y_range;
00054   max_range = ( z_range > max_range ) ? z_range : max_range;
00055   
00056   bboxCenter.x(( min_corner.x() + max_corner.x() ) / 2.0);
00057   bboxCenter.y(( min_corner.y() + max_corner.y() ) / 2.0);
00058   bboxCenter.z(( min_corner.z() + max_corner.z() ) / 2.0);
00059   
00060   epsilonBBox = max_range / pow(2.0, OCTREE_MAX_DEPTH_OCTREE[OCTREE_TIME_ACCURACY_LEVEL] ) * 0.15;
00061   bboxDimension = max_range + 2 * epsilonBBox;
00062   
00063   
00064   bboxMaxX = bboxCenter.x() + bboxDimension / 2.0;
00065   bboxMinX = bboxCenter.x() - bboxDimension / 2.0;
00066   bboxMaxY = bboxCenter.y() + bboxDimension / 2.0;
00067   bboxMinY = bboxCenter.y() - bboxDimension / 2.0;
00068   bboxMaxZ = bboxCenter.z() + bboxDimension / 2.0;
00069   bboxMinZ = bboxCenter.z() - bboxDimension / 2.0;
00070   
00071     // Data structure flags are used to know if a perticualr data structure is already established or not. 
00072 
00073   statusCubitOctreeColoring = CUBIT_FALSE;
00074   
00075   minDepth = OCTREE_MIN_DEPTH_OCTREE[OCTREE_TIME_ACCURACY_LEVEL]; 
00076   maxDepth = OCTREE_MAX_DEPTH_OCTREE[OCTREE_TIME_ACCURACY_LEVEL]; 
00077   
00078   root = NULL;
00079 }
00080 
00081 CubitOctree::~CubitOctree(){
00082   
00083   greyCellList.clean_out();
00084   
00085   boundaryWhiteGridNodeList.clean_out();
00086 
00087   int i;
00088   int num_elm;
00089   
00090     //PRINT_INFO("BEFORE: Deleted grid nodes \n");
00091     //system("pause");
00092   num_elm = gridNodeVector.size();
00093   CubitOctreeNode *ognode;
00094   for( i = 0; i < num_elm; i++ ){
00095     ognode = gridNodeVector.pop();
00096     if (ognode) {delete ognode;}
00097   }
00098   for (i=0; i < greyGridNodeVector.size(); ++i)
00099   {
00100     delete (greyGridNodeVector.get_and_step());
00101   }
00102   
00103     //PRINT_INFO("AFTER: Deleted grid nodes \n");
00104     //system("pause");
00105   delete root;
00106 
00107     //PRINT_INFO("AFTER: Deleted CubitOctree cells \n");
00108     // system("pause");
00109 }
00110 
00111 
00112 double CubitOctree::calculate_depth_based_on_size( double size ){
00113 
00114   return ceil(  log((double)get_bbox_dimension() / (double)size ) / log(2.0)  );
00115 }
00116 
00117 
00118 
00119 CubitBoolean CubitOctree::initialize_octree_generation( void ){
00120   
00121   
00122     // Normalize the center and size of the CubitOctree cell
00123     // bbox dimension is 1.0 with center at 1/2, 1/2, 1/2
00124   CubitVector norm_center;
00125   int level = 0;
00126   int i,j, k, l, m, n;
00127   double norm_dim;
00128   
00129     //norm_center.x(100.0 / 2);
00130     //norm_center.y(100.0 / 2);
00131     //norm_center.z(100.0 / 2);
00132     //norm_dim = 100.0;
00133   
00134   norm_center = bboxCenter;
00135   norm_dim = bboxDimension;
00136   
00137     //VERBOSE PRINT_INFO(" BBOX Center = [ %f %f %f ] Dim = %f \n",norm_center.x(), norm_center.y(), norm_center.z(), norm_dim );
00138   
00139   root = new CubitOctreeCell( norm_center, norm_dim, level, NULL );
00140   CubitVector shift;
00141   
00142   
00143   
00144   CubitOctreeNode *ptr_oct_node[2][2][2];
00145   
00146   for( i = 0; i < 2; i++ ){
00147     for( j = 0; j < 2; j++ ){
00148       for( k = 0; k < 2; k++ ){
00149         if( i == 0 )
00150           shift.x( -norm_dim / 2.0 );
00151         else
00152           shift.x( norm_dim / 2.0 );
00153           
00154         if( j == 0 )
00155           shift.y( -norm_dim / 2.0 );
00156         else
00157           shift.y( norm_dim / 2.0 );        
00158           
00159         if( k == 0 )
00160           shift.z( -norm_dim / 2.0 );
00161         else
00162           shift.z( norm_dim / 2.0 );
00163           
00164         l = ( i == 1 )?0:1;  // fliped 0 -> 1;  1 -> 0
00165         m = ( j == 1 )?0:1;
00166         n = ( k == 1 )?0:1;
00167           
00168         ptr_oct_node[i][j][k] = new CubitOctreeNode( norm_center + shift, root, l, m, n ); 
00169         apd_vector_item( ptr_oct_node[i][j][k] );
00170       }
00171     }
00172   }
00173   
00174     // Connecting CubitOctree Grid Nodes
00175     //- adj_node[0] = O_FRONT Node
00176     //  adj_node[1] = O_BACK Node
00177     //  adj_node[2] = O_RIGHT Node
00178     //  adj_node[3] = O_LEFT Node
00179     //  adj_node[4] = O_TOP Node
00180     //  adj_node[5] = O_BOTTOM Node
00181   
00182   for( i = 0; i < 2; i++ ){
00183     for( j = 0; j < 2; j++ ){
00184       for( k = 0; k < 2; k++ ){
00185     
00186         if( i == 0 ){
00187           ptr_oct_node[i][j][k]->set_adj_node( O_FRONT, ptr_oct_node[1][j][k] );
00188           ptr_oct_node[i][j][k]->set_adj_node_distance( O_FRONT, 0 );
00189         }
00190         else{
00191           ptr_oct_node[i][j][k]->set_adj_node( O_BACK, ptr_oct_node[0][j][k] );
00192           ptr_oct_node[i][j][k]->set_adj_node_distance( O_BACK, 0 );
00193         }
00194           
00195         if( j == 0 ){
00196           ptr_oct_node[i][j][k]->set_adj_node( O_RIGHT, ptr_oct_node[i][1][k] );
00197           ptr_oct_node[i][j][k]->set_adj_node_distance( O_RIGHT, 0 );
00198         }
00199         else{
00200           ptr_oct_node[i][j][k]->set_adj_node( O_LEFT, ptr_oct_node[i][0][k] );
00201           ptr_oct_node[i][j][k]->set_adj_node_distance( O_LEFT, 0 );
00202         }
00203           
00204         if( k == 0 ){
00205           ptr_oct_node[i][j][k]->set_adj_node( O_TOP, ptr_oct_node[i][j][1] );
00206           ptr_oct_node[i][j][k]->set_adj_node_distance( O_TOP, 0 );
00207         }
00208         else{
00209           ptr_oct_node[i][j][k]->set_adj_node( O_BOTTOM, ptr_oct_node[i][j][0] );
00210           ptr_oct_node[i][j][k]->set_adj_node_distance( O_BOTTOM, 0 );
00211         }
00212     
00213       }
00214     }
00215   }
00216   
00217   
00218     // Update the nodes of the root
00219   for( i = 0; i < 2; i++ ){
00220     for( j = 0; j < 2; j++ ){
00221       for( k = 0; k < 2; k++ ){
00222         root->set_oct_grid_node( i, j, k, ptr_oct_node[i][j][k] );
00223     
00224       }
00225     }
00226   }
00227   
00228     //root->display( this );
00229   
00230   return CUBIT_TRUE;
00231 }
00232 
00233 
00234 CubitBoolean CubitOctree::build_octree_till_min_depth( CubitOctreeCell *ptr_cell ){
00235   
00236   int i, j, k;
00237   
00238   if( ptr_cell->get_depth() < minDepth ){
00239     
00240     subdivide_octree_cell( ptr_cell, NULL, minDepth );
00241     
00242       // PRINT_INFO(" establish connection between new nodes and pointed old nodes \n");
00243     
00244       // Recursively call all the childen       
00245     for( i = 0; i < 2; i++ ){
00246       for( j = 0; j < 2; j++ ){
00247         for( k = 0; k < 2; k++ ){
00248             // PRINT_INFO("Child : %d %d %d \n",i,j,k );
00249           build_octree_till_min_depth( ptr_cell->get_child(i, j, k) );
00250         }
00251       }
00252     }
00253 
00254       // PRINT_INFO(" finished calling all recursively all children \n");
00255   
00256   } 
00257   return CUBIT_TRUE;
00258   
00259 }
00260 
00261 
00262 
00263 CubitBoolean CubitOctree::subdivide_octree_based_on_facet_point( OctreeFacetPointData *ptr_facet_point_data, int max_depth  ){
00264   
00265   CubitOctreeCell *ptr_leaf_cell;
00266 
00267   ptr_leaf_cell = root->find_leaf_octree_cell( ptr_facet_point_data->coordinates() );
00268   
00269   CubitBoolean result = subdivide_octree_from_leaf_cell( ptr_leaf_cell, ptr_facet_point_data, max_depth );
00270   if( result == CUBIT_FALSE ){
00271       // Simillar point exist; delete this facet_point_data
00272     delete ptr_facet_point_data;
00273     return CUBIT_FALSE;
00274   }
00275 
00276   return CUBIT_TRUE;
00277 }
00278 
00279 
00280 CubitBoolean CubitOctree::find_cells_based_on_surface_angle(DLIList<CubitOctreeCell*> &refine_cell_list, DLIList<CubitFacetEdge*> &crease_edges, DLIList<CubitFacet*> &crease_facets,
00281                                                        RefFace *one, RefFace *two, const CubitBoolean draw_facets, double dihedral_angle_thresh)
00282 {
00283   int i, j, k;
00284   
00285   CubitFacetEdge *edge;
00286   CubitBoolean water_tight = CUBIT_FALSE;
00287   
00288   DLIList<DLIList<CubitFacet*>*> shell_list;
00289   DLIList<CubitQuadFacet*> dummy;
00290   DLIList<CubitFacet*> mark1, mark2, facets, facet1, facet2, draw_facet_list;
00291   DLIList<CubitFacetEdge*> new_edge1, new_edge2;
00292   DLIList<CubitPoint*> new_points;
00293 
00294   CubitFacet *temp_facet = NULL;
00295   
00296   if (one == NULL || two == NULL)
00297   {
00298     return CUBIT_FALSE;
00299   }
00300 
00301   dihedral_angle_thresh = cos(dihedral_angle_thresh);
00302 
00303   TDOctreeRefFace *td_skl_one = TDOctreeRefFace::get_td(one),
00304       *td_skl_two = TDOctreeRefFace::get_td(two);
00305 
00306   if (td_skl_one == NULL || td_skl_two == NULL ||
00307       td_skl_one->get_ptr_cubit_facet_list() == NULL || td_skl_two->get_ptr_cubit_facet_list() == NULL)
00308   {
00309     PRINT_DEBUG_157("CubitOctree::find_cells_based_on_surface_angle: null TDs or null facet list ptr\n");
00310     PRINT_DEBUG_157("RefFace one->id()=%d, td_skl=%p\ntwo->id()=%d, td_skl=%p\n", one->id(), (void*)td_skl_one,
00311                     two->id(), (void*)td_skl_two);
00312     PRINT_DEBUG_157("This may happen when one volume sizing function \'steals\'\na surface from an existing sizing function.");
00313     return CUBIT_FALSE;
00314   }
00315 
00316   if (!td_skl_one->get_create_2dmat() || !td_skl_two->get_create_2dmat()) {return CUBIT_FALSE;}
00317   
00318 
00319   FacetDataUtil::copy_facets(*TDOctreeRefFace::get_td(one)->get_ptr_cubit_facet_list(), facet1, new_points, new_edge1);
00320   FacetDataUtil::copy_facets(*TDOctreeRefFace::get_td(two)->get_ptr_cubit_facet_list(), facet2, new_points, new_edge2);
00321 
00322   for (i=0; i < new_edge1.size(); ++i)
00323   {
00324     edge = new_edge1.get_and_step();
00325     if (edge->number_tris() == 1)
00326     {
00327       temp_facet = edge->adj_facet(0);
00328       if (!temp_facet->marked())
00329       {
00330         mark1.append(edge->adj_facet(0));
00331         temp_facet->marked(CUBIT_TRUE);
00332       }
00333         //facets.append(edge->adj_facet(0));
00334     }
00335   }
00336 
00337   for (i=0; i < new_edge2.size(); ++i)
00338   {
00339     edge = new_edge2.get_and_step();
00340     if (edge->number_tris() == 1)
00341     {
00342       temp_facet = edge->adj_facet(0);
00343       if (!temp_facet->marked())
00344       {
00345         mark2.append(edge->adj_facet(0));
00346         temp_facet->marked(CUBIT_TRUE);
00347       }
00348         //facets.append(edge->adj_facet(0));
00349     }
00350   }
00351 
00352   for (i=0; i < mark1.size(); ++i)
00353   {
00354     mark1.get_and_step()->marked(CUBIT_FALSE);
00355   }
00356   
00357   for (i=0; i < mark2.size(); ++i)
00358   {
00359     mark2.get_and_step()->marked(CUBIT_FALSE);
00360   }
00361   
00362   facets += facet1;
00363   facets += facet2;
00364   
00365   facet1.clean_out(); //delete facet1;
00366   facet2.clean_out(); //delete facet2;
00367   new_edge1.clean_out(); //delete new_edge1;
00368   new_edge2.clean_out(); //delete new_edge2;
00369   new_points.clean_out(); //delete new_points;
00370   
00371 
00372   CubitStatus rv = FacetDataUtil::split_into_shells(/*facet1*/facets, dummy, shell_list, water_tight);
00373   if (!rv || shell_list.size() == 0)
00374   {
00375     PRINT_DEBUG_157("CubitOctree::find_cells_based_on_surface_angle: could not form shells from facet list!\n");
00376     PRINT_DEBUG_157("return value was %d, number of shells is %d\n", rv, shell_list.size());
00377       //SVDrawTool::clear_non_retained();
00378       //SVDrawTool::draw_facets(facets, CUBIT_RED_INDEX);
00379       //SVDrawTool::mouse_xforms();
00380     return CUBIT_FALSE;
00381   }
00382 
00383       
00384     /*
00385   int num_null_shells = 0;
00386   for (i=0; i < shell_list.size(); ++i)
00387   {
00388     if (shell_list.get_and_step() == NULL)
00389     {
00390       ++num_null_shells;
00391     }
00392   }
00393   if (num_null_shells)
00394   {
00395     PRINT_INFO("CubitOctree::find_cells_based_on_surface_angle: %d of %d shells formed from facets are NULL!\n", num_null_shells, shell_list.size());
00396   }
00397     */
00398       
00400   rv = FacetDataUtil::stitch_facets(shell_list, 1e-4, water_tight, CUBIT_FALSE);
00401   if (!rv)
00402   {
00403     PRINT_DEBUG_157("CubitOctree::find_cells_based_on_surface_angle: could not stitch facets!\n");
00404       //SVDrawTool::clear_non_retained();
00405       //SVDrawTool::draw_facets(facets, CUBIT_RED_INDEX);
00406       //SVDrawTool::mouse_xforms();
00407     return CUBIT_FALSE;
00408   }
00410   
00411     //TDSkeletonCubitFacet::stitch_facets(facet1,facet2);
00412       
00413   for (i=0; i < mark1.size(); ++i)
00414   {
00415     CubitFacet *facet = mark1.get_and_step();
00416     
00417     if (facet->marked()) {continue;}
00418     facet->marked(CUBIT_TRUE);
00419     
00420     for (j=0; j < mark2.size(); ++j)
00421     {
00422       for (k=0; k < 3; ++k)
00423       {
00424         CubitFacet *other = facet->shared_facet(facet->point(k), facet->point((k+1)%3));
00425         if (other != NULL && other == mark2[j] && !mark2[j]->marked() && (facet->normal() % other->normal() < dihedral_angle_thresh))
00426         {
00427             //crease_edges.append(facet->point(k)->shared_edge(facet->point((k+1)%3)));
00428           find_cell_list_intersecting_line_segment(facet->point(k)->coordinates(), facet->point((k+1)%3)->coordinates(), refine_cell_list);
00429           mark2[j]->marked(CUBIT_TRUE);
00430           
00431           if (draw_facets)
00432           {
00433             draw_facet_list.append(facet);
00434             draw_facet_list.append(other);
00435           }
00436         }
00437       }  
00438     }
00439   }
00440   
00441   if (draw_facets) {
00442     //SVDrawTool::draw_facets(draw_facet_list, CUBIT_MAGENTA_INDEX);
00443   }
00444 
00445   if (refine_cell_list.size() > 0)
00446   {
00447       //crease_facets += facets;
00448     FacetDataUtil::delete_facets(facets);
00449     return CUBIT_TRUE;
00450   }
00451 
00452   FacetDataUtil::delete_facets(facets);
00453   return CUBIT_FALSE;
00454 }
00455 
00456 void CubitOctree::refine_cells_to_target_depth(DLIList<CubitOctreeCell*> &refine_cell_list, const int target_depth)
00457 {
00458   int i, j, k;
00459   CubitOctreeCell *current_cell = NULL;
00460   
00461   while (refine_cell_list.size() > 0)
00462   {
00463     current_cell = refine_cell_list.get();
00464     if (subdivide_octree_cell(current_cell, NULL, target_depth))
00465     {
00466       for (i=0; i < 2; ++i)
00467       {
00468         for (j=0; j < 2; ++j)
00469         {
00470           for (k=0; k < 2; ++k)
00471           {
00472             if (current_cell->get_child(i,j,k) != NULL)
00473             {
00474               refine_cell_list.append(current_cell->get_child(i,j,k));
00475             }
00476           }
00477         }
00478       }
00479     }
00480     refine_cell_list.extract();
00481   }
00482 }
00483 
00484 // Checks for only marked cells and unmarks it if it doesn't interesect
00485 CubitBoolean CubitOctree::unmark_octree_cells_not_intersecting_with_facet( CubitFacet *ptr_facet, DLIList<CubitOctreeCell *> &CubitOctree_cell_list){
00486 
00487   int i;
00488   CubitOctreeCell *ptr_cell;
00489 
00490   for( i = 0; i < CubitOctree_cell_list.size(); i++ ){
00491     ptr_cell = CubitOctree_cell_list.get_and_step();
00492     if( ptr_cell->get_mark() ==  CUBIT_TRUE ){
00493       if( ptr_cell->does_facet_intersect_octreecell( ptr_facet ) == CUBIT_FALSE ){
00494         ptr_cell->set_mark( CUBIT_FALSE );
00495       }
00496     }
00497   }
00498   return CUBIT_TRUE;
00499 }
00500 
00501 
00502 CubitBoolean CubitOctree::mark_cells_that_intersect_facet_plane( /*CubitFacet *ptr_facet,*/ DLIList<CubitOctreeCell*> &CubitOctree_cell_list )
00503 {
00504   int i;
00505   CubitOctreeCell *ptr_cell;
00506   
00507   for( i = 0; i < CubitOctree_cell_list.size(); i++ ){
00508     ptr_cell = CubitOctree_cell_list.get_and_step();
00509     if( ptr_cell->does_contain_positive_and_negative_nodes() == CUBIT_TRUE ){
00510       ptr_cell->set_mark( CUBIT_TRUE );
00511     }
00512   }
00513    
00514   return CUBIT_TRUE;
00515 }
00516 
00517 
00518 CubitBoolean CubitOctree::subdivide_octree_from_leaf_cell( CubitOctreeCell *ptr_cell, OctreeFacetPointData *facet_point_data, int max_depth ){
00519   
00520   CubitVector cell_center = ptr_cell->get_center();
00521   int i, j, k;
00522   
00523   CubitVector coord;
00524 
00525     // Check if any other coord already present in the cell is close to new facet_point_data
00526     // At the common curves the facet points of adjacent surfaces can come very close to one another
00527     // Therefore id will not work
00528 
00529   if( ptr_cell->num_of_facet_point_data() != 0 ){
00530     if( ptr_cell->is_facet_point_data_present( facet_point_data ) ){
00531       return CUBIT_FALSE;
00532     }
00533   }
00534   
00535 
00536     // If CubitOctree cell contains on fact points just insert the point don't subdivide (Ref. PR-CubitOctree property)
00537     // 
00538   if( ptr_cell->num_of_facet_point_data() == 0 || ptr_cell->get_depth() >= max_depth ){
00539       
00540     ptr_cell->append_list_item( facet_point_data );
00541     return CUBIT_TRUE;
00542     
00543   }
00544 
00545   coord = facet_point_data->coordinates();  
00546   bool outside_tolerance = CUBIT_TRUE;
00547   DLIList<OctreeFacetPointData*> *fpdata = ptr_cell->get_facet_point_data_list();
00548   for (i=0; i < fpdata->size(); ++i)
00549   {
00550     if ((fpdata->get_and_step()->coordinates()-coord).length_squared() <= OCTREE_EPSILON * OCTREE_EPSILON)
00551     {
00552       outside_tolerance = CUBIT_FALSE;
00553     }
00554   }
00555   
00556   if (!outside_tolerance)
00557   {
00558     ptr_cell->append_list_item(facet_point_data);
00559     return CUBIT_TRUE;
00560     
00561   }
00562   else{
00563   
00564     subdivide_octree_cell( ptr_cell, NULL, max_depth );
00565     
00566     
00567       // find the leaf cell containing the point
00568       // call the subdivision recursively from leaf cell containing the point.
00569     i = j = k = 1;
00570     
00571     if( coord.x() < cell_center.x() )
00572       i = 0; 
00573     
00574     if( coord.y() < cell_center.y() )
00575       j = 0;
00576     
00577     if( coord.z() < cell_center.z() )
00578       k = 0; 
00579     
00580     subdivide_octree_from_leaf_cell( ptr_cell->get_child(i, j, k), facet_point_data, max_depth );
00581     
00582   }
00583   return CUBIT_TRUE;
00584   
00585 }
00586 
00587 CubitBoolean CubitOctree::find_octree_cells_contained_inside_bbox( CubitFacet *ptr_facet, DLIList<CubitOctreeCell*> &CubitOctree_cell_list ){
00588   
00589   CubitPoint *p0, *p1, *p2;
00590   CubitBox facet_bbox;
00591   DLIList<CubitOctreeCell*> queue;
00592   DLIList<CubitOctreeCell*> marked_list;
00593   CubitOctreeCell *CubitOctree_cell1, *CubitOctree_cell2, *CubitOctree_cell3;
00594   CubitOctreeCell *ptr_cell;
00595   CubitBoolean result;
00596   int i;
00597   
00598   queue.clean_out();
00599   
00600   facet_bbox = ptr_facet->bounding_box();
00601     //PRINT_DEBUG_157(" BBox min : %f %f %f max : %f %f %f ", facet_bbox.minimum().x(), facet_bbox.minimum().y(), facet_bbox.minimum().z(), facet_bbox.maximum().x(),  facet_bbox.maximum().y(),  facet_bbox.maximum().z() );
00602   
00603   ptr_facet->points( p0, p1, p2 );
00604   
00605   
00606   CubitOctree_cell1 = find_octreecell( p0->coordinates() );
00607     //CubitOctree_cell_list.push( CubitOctree_cell1 );
00608   queue.push( CubitOctree_cell1 );
00609   CubitOctree_cell1->set_mark( CUBIT_TRUE );
00610   
00611   CubitOctree_cell2 = find_octreecell( p1->coordinates() );
00612     //CubitOctree_cell_list.push( CubitOctree_cell2 );
00613   if( CubitOctree_cell2->get_mark() == CUBIT_FALSE ){
00614     queue.push( CubitOctree_cell2 );
00615     CubitOctree_cell2->set_mark( CUBIT_TRUE );
00616   }
00617   
00618   CubitOctree_cell3 = find_octreecell( p2->coordinates() );
00619     //CubitOctree_cell_list.push( CubitOctree_cell3 );
00620   if( CubitOctree_cell3->get_mark() == CUBIT_FALSE ){
00621     queue.push( CubitOctree_cell3 );
00622     CubitOctree_cell3->set_mark( CUBIT_TRUE );
00623   }
00624   
00625   while( queue.size() > 0 ){
00626     
00627     ptr_cell = queue.pop();
00628     marked_list.push( ptr_cell );
00629     
00630     result = ptr_cell->is_intersects_box( facet_bbox );     
00631     
00632     if( result == CUBIT_TRUE ){
00633       CubitOctree_cell_list.push( ptr_cell );
00634         // add adjacent unmarked CubitOctree cells to queue
00635       ptr_cell->add_adjacent_unmarked_cells( queue );
00636     }
00637     
00638   }
00639   
00640   for( i = 0; i < marked_list.size(); i++ ){        
00641     ptr_cell = marked_list.get_and_step();
00642     ptr_cell->set_mark( CUBIT_FALSE );
00643       //ptr_cell->display(this);
00644   }
00645   
00646   return CUBIT_TRUE;
00647 }
00648 
00649 
00650 CubitBoolean CubitOctree::mark_positive_and_negative_octree_grid_nodes( CubitFacet *ptr_facet, DLIList<CubitOctreeCell*> &CubitOctree_cell_list, DLIList<CubitOctreeNode *> &CubitOctree_grid_node_list ){
00651   
00652   CubitOctreeCell *ptr_cell;
00653   int i, l, m, n;
00654   CubitOctreeNode *ptr_grid_node;
00655   
00656   for( i = 0; i < CubitOctree_cell_list.size(); i++ ){
00657     ptr_cell = CubitOctree_cell_list.get_and_step();
00658     
00659     for( l = 0; l < 2; l++ ){
00660       for( m = 0; m < 2; m++ ){
00661         for( n = 0; n < 2; n++ ){
00662           ptr_grid_node = ptr_cell->get_octree_grid_node( l, m, n );
00663       
00664           if( ptr_grid_node->get_mark() == CUBIT_FALSE ){
00665             ptr_grid_node->find_half_space( ptr_facet );
00666             ptr_grid_node->set_mark( CUBIT_TRUE );
00667             CubitOctree_grid_node_list.push( ptr_grid_node );
00668           }
00669         }
00670       }
00671     }
00672   }
00673   
00674   return CUBIT_TRUE;
00675 }
00676 
00677 
00678 CubitBoolean CubitOctree::find_intersection_between_grid_edges_and_facet(  CubitOctreeType type, RefFace *ptr_face, CubitFacet *ptr_facet, DLIList<CubitOctreeNode *> &CubitOctree_grid_node_list  ){
00679   
00680   int i;
00681   CubitOctreeNode *ptr_grid_node;
00682   
00683   for( i = 0; i < CubitOctree_grid_node_list.size(); i++ ){
00684     ptr_grid_node = CubitOctree_grid_node_list.get_and_step();
00685     if( ptr_grid_node->get_halfspace_direction() == OCTREE_POSITIVE ){
00686       ptr_grid_node->find_intersection_with_facet( type,  ptr_face, ptr_facet, boundaryWhiteGridNodeList );
00687     }
00688   }
00689   
00690     // reset all +ve grid  nodes to -ve
00691   for( i = 0; i < CubitOctree_grid_node_list.size(); i++ ){
00692     ptr_grid_node = CubitOctree_grid_node_list.get_and_step();
00693     ptr_grid_node->set_mark( CUBIT_FALSE ); 
00694     if( ptr_grid_node->get_halfspace_direction() == OCTREE_POSITIVE ){
00695       ptr_grid_node->set_halfspace_direction( OCTREE_NEGATIVE );
00696     }
00697   }
00698   
00699   return CUBIT_TRUE;
00700 }
00701 
00702 
00703 void CubitOctree::gather_initial_grid_nodes_for_smooth_transition( DLIList<CubitOctreeNode*> &queue ){
00704   int i; 
00705   CubitOctreeNode *ptr_grid_node;
00706   int depth_difference;
00707 
00708   for( i = 0; i < gridNodeVector.size(); i++ ){
00709     ptr_grid_node = gridNodeVector.get_and_step();
00710     
00711     depth_difference = ptr_grid_node->find_min_depth_cell_and_depth_difference();
00712       // By default mark is maintained CUBIT_FALSE
00713     if( depth_difference >= 2 ){
00714       queue.push( ptr_grid_node );
00715       ptr_grid_node->set_mark( CUBIT_TRUE );  
00716     }
00717   }
00718 
00719 }
00720 
00721 
00722 // Visit every black node and check the depth of the adjacent cells.
00723 // If the cells differ more than 1 then split the maximum cell.
00724 CubitBoolean CubitOctree::establish_smooth_transition_of_cells( int max_depth ){
00725   
00726   CubitOctreeNode *ptr_grid_node;
00727   int depth_difference;
00728   DLIList<CubitOctreeNode *> queue;
00729 
00730   gather_initial_grid_nodes_for_smooth_transition( queue );
00731 
00732   while( queue.size() > 0 ){
00733     ptr_grid_node = queue.remove();
00734     ptr_grid_node->set_mark( CUBIT_FALSE );
00735       // WARNING: Why Cubit Black only
00736       // if( ptr_grid_node->get_color() == CUBIT_BLACK_INDEX ){
00737     depth_difference = ptr_grid_node->get_cell_depth_difference();
00738     if( depth_difference >= 2 ){
00739       subdivide_octree_cell( ptr_grid_node->get_min_depth_cell(), &queue, max_depth );
00740     }
00741       //}
00742   }
00743 
00744   return CUBIT_TRUE;
00745 }
00746 
00747 CubitBoolean CubitOctree::subdivide_octree_cell( CubitOctreeCell *ptr_cell, DLIList<CubitOctreeNode*> *ptr_queue, int max_depth ){
00748 
00749   CubitVector cell_center, shift;
00750   int new_depth, i, j, k, l, m, n, p, q, r;
00751   double norm_dim, total_dist, new_half_point;
00752   CubitOctreeCell *ptr_child_cell;
00753   CubitOctreeNode *local_node_matrix[3][3][3], *ptr_node, *ptr_new_oct_grid_node;
00754   CubitBoolean new_node_matrix[3][3][3];
00755   CubitOctreeCell *ptr_new_cell;
00756 
00757   if( ptr_cell->get_depth() < max_depth && ptr_cell->is_leaf() )
00758   {
00759     ptr_cell->set_leaf( CUBIT_FALSE );
00760     norm_dim = ptr_cell->get_dimension();
00761     new_depth = ptr_cell->get_depth() + 1;
00762     cell_center = ptr_cell->get_center();
00763       // PRINT_INFO("Num = %d Depth = %d  Coord = %2.12e %2.12e %2.12e\n", ptr_cell->get_num(), depth, cell_center.x(), cell_center.y(), cell_center.z() );
00764         
00765     for( i = 0; i < 2; i++ ){
00766       for( j = 0; j < 2; j++ ){
00767         for( k = 0; k < 2; k++ ){
00768             
00769           if( i == 0 )
00770             shift.x( -norm_dim / 4.0 );
00771           else
00772             shift.x( norm_dim / 4.0 );
00773             
00774           if( j == 0 )
00775             shift.y( -norm_dim / 4.0 );
00776           else
00777             shift.y( norm_dim / 4.0 );      
00778             
00779           if( k == 0 )
00780             shift.z( -norm_dim / 4.0 );
00781           else
00782             shift.z( norm_dim / 4.0 );
00783             
00784             
00785           ptr_child_cell = new CubitOctreeCell( cell_center  + shift, norm_dim / 2.0, new_depth, ptr_cell );
00786           ptr_cell->set_child( i, j, k, ptr_child_cell );
00787             
00788         }
00789       }
00790     }
00791     
00792       // PRINT_INFO(" 8 children created \n");
00793       // Update the wire frame model 
00794     
00795     for( i = 0; i < 3; i++ ){
00796       for( j = 0; j < 3; j++ ){
00797         for( k = 0; k < 3; k++ ){
00798           new_node_matrix[i][j][k] = CUBIT_FALSE;
00799         }
00800       }
00801     }
00802     
00803       // find the eight corner nodes of the local matrix
00804     for( i = 0; i < 2; i++ ){
00805       for( j = 0; j < 2; j++ ){
00806         for( k = 0; k < 2; k++ ){
00807           if( i == 1 )
00808             l = 2;
00809           else
00810             l = 0;
00811             
00812           if( j == 1 )
00813             m = 2;
00814           else
00815             m = 0;
00816             
00817           if( k == 1 )
00818             n = 2;
00819           else
00820             n = 0;
00821             
00822           local_node_matrix[l][m][n] = ptr_cell->get_octree_grid_node( i, j, k );
00823         }
00824       }
00825     }
00826     
00827       // PRINT_INFO(" find 8 corner nodes of local 3 x 3 matrix \n");
00828 
00829     new_half_point = 1.0 / pow(2.0, new_depth);
00830     
00831       // To find 8
00832     ptr_node = local_node_matrix[2][0][2];
00833     total_dist = 0.0;
00834     while(  total_dist < new_half_point ){
00835       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BOTTOM ) );
00836       ptr_node = ptr_node->get_adj_node( O_BOTTOM );            
00837     }
00838     
00839     if( total_dist == new_half_point ){
00840       local_node_matrix[2][0][1] = ptr_node;
00841     }
00842     else{
00843       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][2])->x(), (local_node_matrix[2][0][2])->y(), (local_node_matrix[2][0][2])->z() - norm_dim /2 );
00844       apd_vector_item( ptr_new_oct_grid_node );
00845       local_node_matrix[2][0][1] = ptr_new_oct_grid_node;
00846       new_node_matrix[2][0][1] = CUBIT_TRUE;
00847     }
00848     
00849     
00850       // To find 9
00851     ptr_node = local_node_matrix[2][2][2];
00852     total_dist = 0.0;
00853     while( total_dist < new_half_point ){
00854       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BOTTOM ) );
00855       ptr_node = ptr_node->get_adj_node( O_BOTTOM );
00856     }
00857     
00858     if( total_dist == new_half_point ){
00859       local_node_matrix[2][2][1] = ptr_node;
00860     }
00861     else{
00862       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][2][2])->x(), (local_node_matrix[2][2][2])->y(), (local_node_matrix[2][2][2])->z() - norm_dim /2 );
00863       apd_vector_item( ptr_new_oct_grid_node );
00864       local_node_matrix[2][2][1] = ptr_new_oct_grid_node;
00865       new_node_matrix[2][2][1] = CUBIT_TRUE;
00866     }
00867     
00868       // To find 10
00869     ptr_node = local_node_matrix[0][2][2];
00870     total_dist = 0.0;
00871     while( total_dist < new_half_point ){
00872       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BOTTOM ) );
00873       ptr_node = ptr_node->get_adj_node( O_BOTTOM );
00874     }
00875     
00876     if( total_dist == new_half_point ){
00877       local_node_matrix[0][2][1] = ptr_node;
00878     }
00879     else{
00880       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[0][2][2])->x(), (local_node_matrix[0][2][2])->y(), (local_node_matrix[0][2][2])->z() - norm_dim /2 );
00881       apd_vector_item( ptr_new_oct_grid_node );
00882       local_node_matrix[0][2][1] = ptr_new_oct_grid_node;
00883       new_node_matrix[0][2][1] = CUBIT_TRUE;
00884     }
00885     
00886       // To find 11
00887     ptr_node = local_node_matrix[0][0][2];
00888     total_dist = 0.0;
00889     while( total_dist < new_half_point ){
00890       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BOTTOM ) );
00891       ptr_node = ptr_node->get_adj_node( O_BOTTOM );
00892     }       
00893     if( total_dist == new_half_point ){
00894       local_node_matrix[0][0][1] = ptr_node;
00895     }
00896     else{
00897       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[0][0][2])->x(), (local_node_matrix[0][0][2])->y(), (local_node_matrix[0][0][2])->z() - norm_dim /2 );
00898       apd_vector_item( ptr_new_oct_grid_node );
00899       local_node_matrix[0][0][1] = ptr_new_oct_grid_node;
00900       new_node_matrix[0][0][1] = CUBIT_TRUE;
00901     }
00902     
00903       // To find 12
00904     ptr_node = local_node_matrix[2][0][2];
00905     total_dist = 0.0;
00906     while( total_dist < new_half_point ){
00907       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
00908       ptr_node = ptr_node->get_adj_node( O_RIGHT );
00909     }       
00910     if( total_dist == new_half_point ){
00911       local_node_matrix[2][1][2] = ptr_node;
00912     }
00913     else{
00914       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][2])->x(), (local_node_matrix[2][0][2])->y() + norm_dim / 2.0, (local_node_matrix[2][0][2])->z() ); 
00915       apd_vector_item( ptr_new_oct_grid_node );
00916       local_node_matrix[2][1][2] = ptr_new_oct_grid_node;
00917       new_node_matrix[2][1][2] = CUBIT_TRUE;
00918     }
00919     
00920     
00921       // To find 13
00922     ptr_node = local_node_matrix[2][0][0];
00923     total_dist = 0.0;
00924     while( total_dist < new_half_point ){
00925       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
00926       ptr_node = ptr_node->get_adj_node( O_RIGHT );
00927     }       
00928     if( total_dist == new_half_point ){
00929       local_node_matrix[2][1][0] = ptr_node;
00930     }
00931     else{
00932       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][0])->x(), (local_node_matrix[2][0][0])->y() + norm_dim / 2.0, (local_node_matrix[2][0][0])->z() ); 
00933       apd_vector_item( ptr_new_oct_grid_node );
00934       local_node_matrix[2][1][0] = ptr_new_oct_grid_node;
00935       new_node_matrix[2][1][0] = CUBIT_TRUE;
00936     }
00937     
00938       // To find 14
00939     ptr_node = local_node_matrix[0][0][0];
00940     total_dist = 0.0;
00941     while( total_dist < new_half_point ){
00942       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
00943       ptr_node = ptr_node->get_adj_node( O_RIGHT );
00944     }       
00945     if( total_dist == new_half_point ){
00946       local_node_matrix[0][1][0] = ptr_node;
00947     }
00948     else{
00949       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[0][0][0])->x(), (local_node_matrix[0][0][0])->y() + norm_dim / 2.0, (local_node_matrix[0][0][0])->z() ); 
00950       apd_vector_item( ptr_new_oct_grid_node );
00951       local_node_matrix[0][1][0] = ptr_new_oct_grid_node;
00952       new_node_matrix[0][1][0] = CUBIT_TRUE;
00953     }
00954     
00955     
00956       // To find 15
00957     ptr_node = local_node_matrix[0][0][2];
00958     total_dist = 0.0;
00959     while( total_dist < new_half_point ){
00960       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
00961       ptr_node = ptr_node->get_adj_node( O_RIGHT );
00962     }       
00963     if( total_dist == new_half_point ){
00964       local_node_matrix[0][1][2] = ptr_node;
00965     }
00966     else{
00967       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[0][0][2])->x(), (local_node_matrix[0][0][2])->y() + norm_dim / 2.0, (local_node_matrix[0][0][2])->z() ); 
00968       apd_vector_item( ptr_new_oct_grid_node );
00969       local_node_matrix[0][1][2] = ptr_new_oct_grid_node;
00970       new_node_matrix[0][1][2] = CUBIT_TRUE;
00971     }
00972     
00973     
00974       // To find 16
00975     ptr_node = local_node_matrix[2][0][2];
00976     total_dist = 0.0;
00977     while( total_dist < new_half_point ){
00978       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
00979       ptr_node = ptr_node->get_adj_node( O_BACK );
00980     }
00981     
00982     if( total_dist == new_half_point ){
00983       local_node_matrix[1][0][2] = ptr_node;
00984     }
00985     else{
00986       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][2])->x() - norm_dim / 2.0, (local_node_matrix[2][0][2])->y(), (local_node_matrix[2][0][2])->z());
00987       apd_vector_item( ptr_new_oct_grid_node );
00988       local_node_matrix[1][0][2] = ptr_new_oct_grid_node;
00989       new_node_matrix[1][0][2] = CUBIT_TRUE;
00990     }
00991     
00992     
00993       // To find 17
00994     ptr_node = local_node_matrix[2][0][0];
00995     total_dist = 0.0;
00996     while( total_dist < new_half_point ){
00997       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
00998       ptr_node = ptr_node->get_adj_node( O_BACK );
00999     }
01000     if( total_dist == new_half_point ){
01001       local_node_matrix[1][0][0] = ptr_node;
01002     }
01003     else{
01004       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][0])->x() - norm_dim / 2.0, (local_node_matrix[2][0][0])->y(), (local_node_matrix[2][0][0])->z());
01005       apd_vector_item( ptr_new_oct_grid_node );
01006       local_node_matrix[1][0][0] = ptr_new_oct_grid_node;
01007       new_node_matrix[1][0][0] = CUBIT_TRUE;
01008     }
01009     
01010     
01011       // To find 18
01012     ptr_node = local_node_matrix[2][2][0];
01013     total_dist = 0.0;
01014     while( total_dist < new_half_point ){
01015       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
01016       ptr_node = ptr_node->get_adj_node( O_BACK );
01017     }
01018     if( total_dist == new_half_point ){
01019       local_node_matrix[1][2][0] = ptr_node;
01020     }
01021     else{
01022       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][2][0])->x() - norm_dim / 2.0, (local_node_matrix[2][2][0])->y(), (local_node_matrix[2][2][0])->z());
01023       apd_vector_item( ptr_new_oct_grid_node );
01024       local_node_matrix[1][2][0] = ptr_new_oct_grid_node;
01025       new_node_matrix[1][2][0] = CUBIT_TRUE;
01026     }
01027     
01028     
01029       // To find 19
01030     ptr_node = local_node_matrix[2][2][2];
01031     total_dist = 0.0;
01032     while( total_dist < new_half_point ){
01033       total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
01034       ptr_node = ptr_node->get_adj_node( O_BACK );
01035     }
01036     if( total_dist == new_half_point ){
01037       local_node_matrix[1][2][2] = ptr_node;
01038     }
01039     else{
01040       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][2][2])->x() - norm_dim / 2.0, (local_node_matrix[2][2][2])->y(), (local_node_matrix[2][2][2])->z());
01041       apd_vector_item( ptr_new_oct_grid_node );
01042       local_node_matrix[1][2][2] = ptr_new_oct_grid_node;
01043       new_node_matrix[1][2][2] = CUBIT_TRUE;
01044     }
01045     
01046       // PRINT_INFO(" found nodes on edges of local 3 x 3 matrix \n");
01047       // To find 20
01048     ptr_node = local_node_matrix[2][0][1];
01049     if( new_node_matrix[2][0][1] == CUBIT_FALSE && ptr_node->get_adj_node( O_RIGHT ) != NULL ){
01050       
01051       total_dist = 0.0;
01052       while( total_dist < new_half_point ){
01053         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
01054         ptr_node = ptr_node->get_adj_node( O_RIGHT );
01055       }
01056       if( total_dist == new_half_point ){
01057         local_node_matrix[2][1][1] = ptr_node;
01058           // PRINT_INFO(" old node 20 found \n");
01059       }
01060       else{
01061         PRINT_INFO("WARNING: Node not found \n");
01062       }
01063       
01064     }
01065     else{
01066       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][1])->x(), (local_node_matrix[2][0][1])->y() + norm_dim / 2.0, (local_node_matrix[2][0][1])->z());
01067       apd_vector_item( ptr_new_oct_grid_node );
01068       local_node_matrix[2][1][1] = ptr_new_oct_grid_node;
01069       new_node_matrix[2][1][1] = CUBIT_TRUE;
01070         // PRINT_INFO(" new node 20 created \n");
01071     }
01072     
01073     
01074     
01075       // To find 21
01076     ptr_node = local_node_matrix[0][0][1];
01077     if( new_node_matrix[0][0][1] == CUBIT_FALSE && ptr_node->get_adj_node( O_RIGHT ) != NULL ){
01078       
01079       total_dist = 0.0;
01080       while( total_dist < new_half_point ){
01081         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
01082         ptr_node = ptr_node->get_adj_node( O_RIGHT );
01083       }
01084       if( total_dist == new_half_point ){
01085         local_node_matrix[0][1][1] = ptr_node;
01086           // PRINT_INFO(" old node 21 found \n");
01087       }
01088       else{
01089         PRINT_INFO("WARNING: Node not found \n");
01090       }
01091     }
01092     else{
01093       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[0][0][1])->x(), (local_node_matrix[0][0][1])->y() + norm_dim / 2.0, (local_node_matrix[0][0][1])->z());
01094       apd_vector_item( ptr_new_oct_grid_node );
01095       local_node_matrix[0][1][1] = ptr_new_oct_grid_node;
01096       new_node_matrix[0][1][1] = CUBIT_TRUE;
01097         // PRINT_INFO(" new node 21 created \n");
01098     }
01099     
01100       // To find 22 
01101     ptr_node = local_node_matrix[2][0][1];
01102     if( new_node_matrix[2][0][1] == CUBIT_FALSE && ptr_node->get_adj_node( O_BACK ) != NULL ){
01103       
01104       total_dist = 0.0;
01105       while( total_dist < new_half_point ){
01106         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
01107         ptr_node = ptr_node->get_adj_node( O_BACK );
01108       }
01109       if( total_dist == new_half_point ){
01110         local_node_matrix[1][0][1] = ptr_node;
01111           // PRINT_INFO(" old node 22 found \n");
01112       }
01113       else{
01114         PRINT_INFO("WARNING: Node not found \n");
01115       }
01116     }
01117     else{
01118       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][0][1])->x() - norm_dim / 2.0, (local_node_matrix[2][0][1])->y(), (local_node_matrix[2][0][1])->z());
01119       apd_vector_item( ptr_new_oct_grid_node );
01120       local_node_matrix[1][0][1] = ptr_new_oct_grid_node;
01121       new_node_matrix[1][0][1] = CUBIT_TRUE;
01122         // PRINT_INFO(" new node 22 created \n");
01123     }
01124     
01125       // To find 23 
01126     ptr_node = local_node_matrix[2][2][1];
01127     if( new_node_matrix[2][2][1] == CUBIT_FALSE && ptr_node->get_adj_node( O_BACK ) != NULL ){
01128       
01129       total_dist = 0.0;
01130       while( total_dist < new_half_point ){
01131         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_BACK ) );
01132         ptr_node = ptr_node->get_adj_node( O_BACK );
01133       }
01134       if( total_dist == new_half_point ){
01135         local_node_matrix[1][2][1] = ptr_node;
01136           // PRINT_INFO(" old node 23 found \n");
01137       }
01138       else{
01139         PRINT_INFO("WARNING: Node not found \n");
01140       }
01141     }
01142     else{
01143       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[2][2][1])->x() - norm_dim / 2.0, (local_node_matrix[2][2][1])->y(), (local_node_matrix[2][2][1])->z());
01144       apd_vector_item( ptr_new_oct_grid_node );
01145       local_node_matrix[1][2][1] = ptr_new_oct_grid_node;
01146       new_node_matrix[1][2][1] = CUBIT_TRUE;
01147         // PRINT_INFO(" new node 23 created \n");
01148     }
01149     
01150       // To find 24
01151     ptr_node = local_node_matrix[1][0][0];
01152     if( new_node_matrix[1][0][0] == CUBIT_FALSE && ptr_node->get_adj_node( O_RIGHT ) != NULL ){
01153       
01154       total_dist = 0.0;
01155       while( total_dist < new_half_point ){
01156         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
01157         ptr_node = ptr_node->get_adj_node( O_RIGHT );
01158       }
01159       if( total_dist == new_half_point ){
01160         local_node_matrix[1][1][0] = ptr_node;
01161           // PRINT_INFO(" old node 24 found \n");
01162       }
01163       else{
01164         PRINT_INFO("WARNING: Node not found \n");
01165       }
01166     }
01167     else{
01168       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[1][0][0])->x(), (local_node_matrix[1][0][0])->y() + norm_dim / 2.0, (local_node_matrix[1][0][0])->z());
01169       apd_vector_item( ptr_new_oct_grid_node );
01170       local_node_matrix[1][1][0] = ptr_new_oct_grid_node;
01171       new_node_matrix[1][1][0] = CUBIT_TRUE;
01172         // PRINT_INFO(" new node 24 created \n");
01173     }
01174     
01175     
01176     
01177       // To find 25
01178     ptr_node = local_node_matrix[1][0][2];
01179     if( new_node_matrix[1][0][2] == CUBIT_FALSE && ptr_node->get_adj_node( O_RIGHT ) != NULL ){
01180       
01181       total_dist = 0.0;
01182       while( total_dist < new_half_point ){
01183         total_dist += 1.0 / pow( 2.0, ptr_node->get_adj_node_distance( O_RIGHT ) );
01184         ptr_node = ptr_node->get_adj_node( O_RIGHT );
01185       }
01186       if( total_dist == new_half_point ){
01187         local_node_matrix[1][1][2] = ptr_node;
01188           // PRINT_INFO(" old node 25 found \n");
01189       }
01190       else{
01191         PRINT_INFO("WARNING: Node not found \n");
01192       }
01193     }
01194     else{
01195       ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[1][0][2])->x(), (local_node_matrix[1][0][2])->y() + norm_dim / 2.0, (local_node_matrix[1][0][2])->z());
01196       apd_vector_item( ptr_new_oct_grid_node );
01197       local_node_matrix[1][1][2] = ptr_new_oct_grid_node;
01198       new_node_matrix[1][1][2] = CUBIT_TRUE;
01199         // PRINT_INFO(" new node 25 created \n");
01200     }
01201     
01202       // PRINT_INFO(" found nodes on faces of local 3 x 3 matrix \n");
01203     
01204       // To find 26
01205     ptr_new_oct_grid_node = new CubitOctreeNode( (local_node_matrix[1][1][2])->x(), (local_node_matrix[1][1][2])->y(), (local_node_matrix[1][1][2])->z() - norm_dim / 2.0);
01206     apd_vector_item( ptr_new_oct_grid_node );
01207     local_node_matrix[1][1][1] = ptr_new_oct_grid_node;
01208     new_node_matrix[1][1][1] = CUBIT_TRUE;
01209     
01210       // PRINT_INFO(" found nodes at center of  local 3 x 3 matrix \n");
01211     
01212       // update the eight nodes of new eight cells
01213     
01214     for( i = 0; i < 2; i++ ){
01215       for( j = 0; j < 2; j++ ){
01216         for( k = 0; k < 2; k++ ){
01217       
01218           ptr_new_cell = ptr_cell->get_child( i, j, k );
01219             
01220           for( l = 0; l < 2; l++ ){
01221             for( m = 0; m < 2; m++ ){
01222               for( n = 0; n < 2; n++ ){
01223                 ptr_new_cell->set_oct_grid_node( l, m, n, local_node_matrix[i+l][j+m][k+n] );
01224     
01225                 p = ( l == 1 )?0:1;  // fliped 0 -> 1;  1 -> 0
01226                 q = ( m == 1 )?0:1;
01227                 r = ( n == 1 )?0:1;
01228                   // update the eight cells of all 27 nodes
01229                 local_node_matrix[i+l][j+m][k+n]->update_adj_cell( ptr_new_cell, p, q, r );
01230               }
01231             }
01232           }
01233         }
01234       }
01235     }
01236 
01237       // PRINT_INFO(" updates nodes of the 8 children and visa versa \n");
01238     
01239       // connnect the new nodes with the old nodes to update the wire frame
01240     
01241     for( i = 0; i < 3; i++ ){
01242       for( j = 0; j < 3; j++ ){
01243         for( k = 0; k < 3; k++ ){
01244           if( new_node_matrix[i][j][k] == CUBIT_TRUE ){
01245             if( i - 1 >= 0 ){
01246               local_node_matrix[i][j][k]->set_adj_node( O_BACK, local_node_matrix[i-1][j][k] );
01247               local_node_matrix[i][j][k]->set_adj_node_distance( O_BACK, new_depth );
01248               local_node_matrix[i-1][j][k]->set_adj_node( O_FRONT, local_node_matrix[i][j][k] );
01249               local_node_matrix[i-1][j][k]->set_adj_node_distance( O_FRONT, new_depth );
01250             }
01251             if( i + 1 <= 2 ){
01252               local_node_matrix[i][j][k]->set_adj_node( O_FRONT, local_node_matrix[i+1][j][k] );
01253               local_node_matrix[i][j][k]->set_adj_node_distance( O_FRONT, new_depth );
01254               local_node_matrix[i+1][j][k]->set_adj_node( O_BACK, local_node_matrix[i][j][k] );             
01255               local_node_matrix[i+1][j][k]->set_adj_node_distance( O_BACK, new_depth );
01256             }
01257               
01258             if( j - 1 >= 0 ){
01259               local_node_matrix[i][j][k]->set_adj_node( O_LEFT, local_node_matrix[i][j-1][k] );
01260               local_node_matrix[i][j][k]->set_adj_node_distance( O_LEFT, new_depth );
01261               local_node_matrix[i][j-1][k]->set_adj_node( O_RIGHT, local_node_matrix[i][j][k] );
01262               local_node_matrix[i][j-1][k]->set_adj_node_distance( O_RIGHT, new_depth );
01263             }
01264             if( j + 1 <= 2 ){
01265               local_node_matrix[i][j][k]->set_adj_node( O_RIGHT, local_node_matrix[i][j+1][k] );
01266               local_node_matrix[i][j][k]->set_adj_node_distance( O_RIGHT, new_depth );
01267               local_node_matrix[i][j+1][k]->set_adj_node( O_LEFT, local_node_matrix[i][j][k] );
01268               local_node_matrix[i][j+1][k]->set_adj_node_distance( O_LEFT, new_depth );
01269             }
01270               
01271             if( k - 1 >= 0 ){
01272               local_node_matrix[i][j][k]->set_adj_node( O_BOTTOM, local_node_matrix[i][j][k-1] );
01273               local_node_matrix[i][j][k]->set_adj_node_distance( O_BOTTOM,  new_depth );
01274               local_node_matrix[i][j][k-1]->set_adj_node( O_TOP, local_node_matrix[i][j][k] );
01275               local_node_matrix[i][j][k-1]->set_adj_node_distance( O_TOP, new_depth );
01276             }
01277             if( k + 1 <= 2 ){
01278               local_node_matrix[i][j][k]->set_adj_node( O_TOP, local_node_matrix[i][j][k+1] );
01279               local_node_matrix[i][j][k]->set_adj_node_distance( O_TOP, new_depth );
01280               local_node_matrix[i][j][k+1]->set_adj_node( O_BOTTOM, local_node_matrix[i][j][k] );
01281               local_node_matrix[i][j][k+1]->set_adj_node_distance( O_BOTTOM, new_depth );
01282             }
01283               
01284           }
01285         }
01286       }
01287     }
01288     
01289       // distribute the facet_points and Ref_face among the children
01290     ptr_cell->distribute_facet_points_among_children();
01291     
01292     
01293     int depth_difference;
01294       // update the depthMaxDepthCell, depthMinDepthCell, cellDepthDifference, and minDepthCell 
01295     if( ptr_queue != NULL ){
01296       for( i = 0; i < 3; i++ ){
01297         for( j = 0; j < 3; j++ ){
01298           for( k = 0; k < 3; k++ ){
01299             depth_difference = local_node_matrix[i][j][k]->find_min_depth_cell_and_depth_difference( );        
01300               /*  RISKY TO OPTIMIZE ANY THING HERE 
01301                   NEEDS CAREFULL STUDY
01302                   if( new_node_matrix[i][j][k] == CUBIT_TRUE ){
01303                     // initialize the variables depthMaxDepthCell, depthMinDepthCell, cellDepthDifference, and minDepthCell  
01304                     local_node_matrix[i][j][k]->find_min_depth_cell_and_depth_difference(  );        
01305                     }
01306                     else{
01307                       // update the variables depthMaxDepthCell, depthMinDepthCell, cellDepthDifference, and minDepthCell 
01308                       max_depth = local_node_matrix[i][j][k]->get_depth_max_depth_cell();
01309                       min_depth = local_node_matrix[i][j][k]->get_depth_min_depth_cell();
01310                       if( new_depth < min_depth || new_depth > max_depth )
01311 
01312                       }
01313               */
01314             if( local_node_matrix[i][j][k]->get_mark() == CUBIT_FALSE  && depth_difference >= 2 ){
01315               ptr_queue->push( local_node_matrix[i][j][k] );
01316               local_node_matrix[i][j][k]->set_mark( CUBIT_TRUE );
01317             }
01318           }
01319         }
01320       }
01321     }
01322 //    ptr_cell->delete_local_node_matrix();
01323     return CUBIT_TRUE;
01324   }
01325   else{
01326     return CUBIT_FALSE;
01327   }
01328 }
01329 
01330 
01331 void CubitOctree::color_octreecell( void ){
01332   
01333   root->coloring(/*greyCellList,*/ blackCellList);
01334   
01335     /* TESTING: Duplicates in grey cell list 
01336        CubitOctreeCell *ptr_cell;
01337        int i;
01338        for( i = 0; i < greyCellList.size(); i++ ){
01339        ptr_cell = greyCellList.get_and_step();
01340        PRINT_INFO("Grey cell ID: %d \n", ptr_cell->id() );
01341        }
01342     */
01343 }
01344 
01345 // p0 and p1 should be within the root's bbox
01346 // note: this is a little sloppy, i'll clean it up when I fix the general neighbor traversal
01347 void CubitOctree::find_cell_list_intersecting_line_segment (const CubitVector &p0, const CubitVector &p1, DLIList<CubitOctreeCell*> &cell_list){
01348 
01349 //  plane order is like in enum OctreePosition
01350   
01351   CubitOctreeCell *start = find_octreecell(p0), *end = find_octreecell(p1);
01352   CubitOctreeCell *current = start, *previous = NULL;
01353   CubitVector curr_pt = p0;
01354   CubitVector dir = p1-p0, center;
01355     //CubitVector ii(1,0,0);
01356     //CubitVector jj(0,1,0);
01357     //CubitVector kk(0,0,1);
01358 
01359   int i;
01360   dir.normalize();
01361   double tx=0, ty=0, tz=0, h=0, a = dir.x(), b = dir.y(), c = dir.z();
01362 
01363   static unsigned int face_corners[6][4] =
01364       {
01365           {0,5,6,2},
01366           {4,3,1,7},
01367           {5,3,7,2},
01368           {0,4,1,6},
01369           {1,7,2,6},
01370           {4,3,5,0}
01371       };
01372 
01373   static unsigned int inwward_cells[6][4][3] =
01374       {
01375           {{0,0,0},{0,0,1},{1,0,0},{1,0,1}},
01376           {{0,1,0},{0,1,1},{1,1,0},{1,1,1}},
01377           {{0,1,0},{0,0,0},{1,0,0},{1,1,0}},
01378           {{0,1,1},{0,0,1},{1,0,1},{1,1,1}},
01379           {{0,0,0},{0,0,1},{0,1,1},{0,1,0}},
01380           {{1,0,0},{1,0,1},{1,1,1},{1,1,0}}
01381       };
01382 
01383   while (current != end)
01384   {
01385     if (current != previous)
01386     {
01387       cell_list.append(current);
01388       previous = current;
01389     }
01390       // do tn, pick face, march stuff here
01391     h = current->get_dimension();
01392     center = current->get_center();
01393     if (fabs(a) > CUBIT_DBL_MIN) {tx = h/(2*fabs(a)) - (curr_pt.x() - center.x())/a;}
01394     else {tx = CUBIT_DBL_MAX;}
01395     if (fabs(b) > CUBIT_DBL_MIN) {ty = h/(2*fabs(b)) - (curr_pt.y() - center.y())/b;}
01396     else {ty = CUBIT_DBL_MAX;}
01397     if (fabs(c) > CUBIT_DBL_MIN) {tz = h/(2*fabs(c)) - (curr_pt.z() - center.z())/c;}
01398     else {tz = CUBIT_DBL_MAX;}
01399 
01400     tx = fabs(tx);
01401     ty = fabs(ty);
01402     tz = fabs(tz);
01403     
01404     double t = 0.0;
01405 
01406     OctreePosition int_plane = O_UNSET;
01407     
01408     if (tx <= ty && tx <= tz)
01409     {
01410       if (a > 0) {t = ((center.x() + h/2) - curr_pt.x())/a; int_plane = O_FRONT;}
01411       else {t = ((center.x() - h/2) - curr_pt.x())/a; int_plane = O_BACK;}
01412     }
01413     
01414     else if (ty <= tx && ty <= tz)
01415     {
01416       if (b > 0) {t = ((center.y() + h/2) - curr_pt.y())/b; int_plane = O_RIGHT;}
01417       else {t = ((center.y() - h/2) - curr_pt.y())/b; int_plane = O_LEFT;}
01418     }
01419     
01420     else if (tz <= ty && tz <= tx)
01421     {
01422       if (c > 0) {t = ((center.z() + h/2) - curr_pt.z())/c; int_plane = O_TOP;}
01423       else {t = ((center.z() - h/2) - curr_pt.z())/c; int_plane = O_BOTTOM;}
01424     }
01425     
01426       // now we have determined the point of intersection on a face of the cell
01427       // based on the face, we should now check the grid nodes' adjacent cells and march if necessary to find the next cell to march in
01428 
01429     CubitOctreeNode* corners[] =
01430         {
01431             current->get_octree_grid_node(1,0,1),
01432             current->get_octree_grid_node(0,1,1),
01433             current->get_octree_grid_node(0,0,0),
01434             current->get_octree_grid_node(1,1,0),
01435           
01436             current->get_octree_grid_node(1,1,1),
01437             current->get_octree_grid_node(1,0,0),
01438             current->get_octree_grid_node(0,0,1),
01439             current->get_octree_grid_node(0,1,0),
01440         };
01441       /*
01442         CubitOctreeNode *cADF = current->get_octree_grid_node(1,0,1); //0
01443         CubitOctreeNode *cBCF = current->get_octree_grid_node(0,1,1); //1
01444         CubitOctreeNode *cCDE = current->get_octree_grid_node(0,0,0); //2
01445         CubitOctreeNode *cABE = current->get_octree_grid_node(1,1,0); //3
01446 
01447         CubitOctreeNode *cABF = current->get_octree_grid_node(1,1,1); //4
01448         CubitOctreeNode *cADE = current->get_octree_grid_node(1,0,0); //5
01449         CubitOctreeNode *cCDF = current->get_octree_grid_node(0,0,1); //6
01450         CubitOctreeNode *cBCE = current->get_octree_grid_node(0,1,0); //7
01451       */
01452       
01453   
01454 
01455       // each digit is an i,j,k index for the inward cell from a grid node
01456 
01457     CubitOctreeNode* int_face[4] = {corners[face_corners[int_plane][0]], corners[face_corners[int_plane][1]], corners[face_corners[int_plane][2]], corners[face_corners[int_plane][3]]};
01458     unsigned int* cells[4] = {inwward_cells[int_plane][0], inwward_cells[int_plane][1], inwward_cells[int_plane][2], inwward_cells[int_plane][3]};
01459 
01460     CubitOctreeCell* four_cells[4];
01461       //CubitOctreeCell *curr_cell = int_face[0]->get_adj_cell(cells[0][0],cells[0][1],cells[0][2]);
01462     bool one_face_neighbor = true;
01463     //bool same_depth = true;
01464     CubitOctreeCell *compare_to = NULL;
01465     
01466     for (i=0; i < 4; ++i)
01467     {
01468       four_cells[i] = int_face[i]->get_adj_cell(cells[i][0],cells[i][1],cells[i][2]);
01469       if (compare_to == NULL && four_cells[i] != NULL) {compare_to = four_cells[i];}
01470     }
01471     
01472     for (i=0; i < 4; ++i)
01473     {
01474       if (four_cells[i] == NULL)
01475       {
01476         one_face_neighbor = false;
01477         //same_depth = false;
01478         break;
01479       }
01480       if (four_cells[i] != NULL)
01481       {
01482         if (four_cells[i] != compare_to) {one_face_neighbor = false;}
01483         //if (four_cells[i]->get_depth() != compare_to->get_depth()) {same_depth = false;}
01484       }
01485     }
01486 
01487     curr_pt += t*dir;
01488     if (one_face_neighbor == true)
01489     {
01490       current = compare_to;//four_cells[0];//curr_cell;.
01491 //      PRINT_INFO("Found one face neighbor for CubitOctree ray traversal\n");
01492     }
01493 
01494     else
01495     {
01496         //curr_pt += (t+.0001)*dir;
01497       int axes[6] = {1,1,2,2,0,0};
01498       int dirs[6] = {-1,1,-1,1,-1,1};
01499       double xyz[3];
01500       curr_pt.get_xyz(xyz);
01501 
01502       xyz[axes[int_plane]] += dirs[int_plane] * (get_bbox_dimension()/(2*pow((double)2.0,maxDepth)));
01503 
01504       CubitVector temp_pt;
01505       temp_pt.set(xyz);
01506         //    SVDrawTool::draw_point(temp_pt, CUBIT_WHITE_INDEX);
01507       
01508       current = find_octreecell(temp_pt);
01509         /* if (current == previous) {
01510            PRINT_INFO("Samet method failed!\n");
01511            }*/
01512       
01513         //PRINT_INFO("Did not find one face neighbor for CubitOctree ray traversal, following Samet method\n");
01514     }
01515       //SVDrawTool::mouse_xforms();
01516     if (current == end) {break;}
01517     
01518   }
01519   cell_list.append(end);
01520   
01521 }
01522 
01523 
01524 
01525 CubitOctreeCell *CubitOctree::find_octreecell( const CubitVector &pnt ){
01526   
01527   return( root->find_leaf_octree_cell( pnt ) );
01528   
01529 }
01530 
01531 CubitBoolean CubitOctree::apd_vector_item( CubitOctreeNode *ptr_node ){
01532   
01533   gridNodeVector.push( ptr_node);
01534   return CUBIT_TRUE;
01535 } 
01536 
01537 
01538 void CubitOctree::find_max_min_size_grid_node_for_scaling( void ){
01539   CubitOctreeNode *ptr_grid_node;
01540   int i;
01541   double size;
01542 
01543   for( i = 0; i < gridNodeVector.size(); i++ ){
01544 
01545     ptr_grid_node = gridNodeVector.get_and_step();
01546     if (ptr_grid_node)
01547     {
01548               
01549       size = ptr_grid_node->get_size(  OCTREE_SIZE_DEFAULT );
01550             
01551       if( size > 0.0 ){
01552           // TESTING
01553           //PRINT_DEBUG_157("size = %f\n", size);
01554               
01555         if( size < minSizeGridNode )
01556           minSizeGridNode = size;
01557               
01558         if( size > maxSizeGridNode )
01559           maxSizeGridNode = size;
01560       }
01561     }
01562             
01563   }
01564 
01565   PRINT_DEBUG_167("Max Size Grid Node =  %f\n", maxSizeGridNode );
01566   PRINT_DEBUG_167("Min Size Grid Node =  %f\n", minSizeGridNode );
01567   double average_size = (minSizeGridNode + maxSizeGridNode) / 2.0;
01568   minSizeGridNode -= .05*average_size;
01569   maxSizeGridNode += .05*average_size;
01570   PRINT_DEBUG_167("For drawing:\n");
01571   PRINT_DEBUG_167("Adjusted Max Size Grid Node =  %f\n", maxSizeGridNode );
01572   PRINT_DEBUG_167("Adjusted Min Size Grid Node =  %f\n", minSizeGridNode );
01573 
01574     //system("pause");
01575     
01576 }
01577 
01578 
01579 
01580 double CubitOctree::size_at_a_point( const CubitVector &point, int size_type ){
01581 //PRINT_DEBUG_116("Entered the SizeAtPoint of SkeletonSizingFunction");
01582   static int mode = 1;
01583   double size = 0.0; 
01584 
01585 //  SVDrawTool::draw_point(point, CUBIT_RED_INDEX);
01586   
01587   
01588   CubitOctreeCell *ptr_cell = find_octreecell( point );
01589   
01590   if( ptr_cell == NULL )
01591   {
01592     
01593     PRINT_INFO("WARNING: Point is outside the solid\n");
01594       //exit(1);
01595     if( mode == 0 )
01596       return size;
01597 
01598       /*
01599         mode = 0;
01600             
01601           // Search the adjacent cells and find the nearest sizing data
01602           CubitVector xyz[6];
01603           CubitVector mod_point;
01604           int i;
01605       
01606           for( i = 0; i < 5; i ++ ){
01607           xyz[i].x(0.0);
01608           xyz[i].y(0.0);
01609           xyz[i].z(0.0);
01610           }
01611       
01612           xyz[0].x( grid_size );
01613           xyz[1].x( -1 * grid_size );
01614           xyz[2].y( grid_size );
01615           xyz[3].y( -1 * grid_size );
01616           xyz[4].z( grid_size );
01617           xyz[5].z( -1 * grid_size );
01618       
01619           double total_size = 0;
01620           int count = 0;
01621           double size1;
01622       
01623           for( i = 0; i < 5; i ++ ){
01624           mod_point = point + xyz[i];
01625           size1 = size_at_point( mod_point, 0.0 );
01626       
01627           if( size1 != 0 ){
01628           total_size += size1;
01629           count++;
01630           }
01631           }
01632       
01633           size = total_size / count;
01634           if( size == 0.0 ){
01635           PRINT_INFO("WARNING: Point is outside the solid can't handle ( size set to average size )\n");
01636           size = 0.5;
01637           size = 2 * ( min_rad + ( size - 1 ) / (10 - 1) * ( max_rad - min_rad ) ) * grid_size * over_all_reduction_scale;
01638           }
01639       
01640           mode = 1;
01641       */
01642   }
01643   else
01644   {
01645     
01646     switch( DEFAULT_INTERPOLATION_INSIDE_CELL ) 
01647     {
01648     case TRILINEAR_INTERPOLATION_INSIDE_CELL:
01649       size = ptr_cell->trilinear_interpolation( point );
01650       break;
01651     
01652     case INVERSE_DISTANCE_INTERPOLATION_INSIDE_CELL:
01653       size = ptr_cell->inverse_distance_interpolation( point );
01654       break;
01655 
01656     case MIN_DISTANCE_INTERPOLATION_INSIDE_CELL:
01657       size = ptr_cell->min_distance_interpolation( point );
01658       break;
01659 
01660     case MIN_SIZE_INTERPOLATION_INSIDE_CELL:
01661       size = ptr_cell->min_size_interpolation( point );
01662       break;
01663 
01664     default:
01665 #ifndef NDEBUG
01666         //PRINT_INFO(" WARNING: Interpolation method inside a cell is not set in SkeletonConstants.hpp\n");
01667 #endif
01668       size = ptr_cell->trilinear_interpolation( point );
01669       break;
01670     }       
01671    
01672     if( size_type == SCALED_SIZE )
01673     {
01674       size = get_scaled_from_wrld_size_grid_node( size );
01675             
01676     }
01677 
01678         
01679       //    if( size_type == MESH_SIZE ){
01680       // Do nothing send size for meshing
01681       //}
01682 
01683   }
01684     // PRINT_INFO("Size at point <%f,%f,%f> = %f\n", point.x(), point.y(), point.z(), size);
01685   return size;
01686 }
01687 
01688 
01689 CubitPointContainment CubitOctree::point_containment( CubitVector point, double tolerance )
01690 {
01691   
01692   CubitOctreeCell *ptr_cell = find_octreecell( point );
01693   
01694   if( ptr_cell == NULL )
01695   {
01696     // not known and needs full acis/mbg point containment
01697     return CUBIT_PNT_BOUNDARY;
01698   }
01699   else
01700   {
01701     switch( ptr_cell->get_color() )
01702     {
01703       case CUBIT_BLACK_INDEX:
01704         return CUBIT_PNT_INSIDE;
01705     
01706       case CUBIT_WHITE_INDEX:
01707         return CUBIT_PNT_OUTSIDE;
01708         
01709       case CUBIT_GREY_INDEX:
01710         return CUBIT_PNT_BOUNDARY;
01711     
01712       default:
01713         return CUBIT_PNT_BOUNDARY;
01714     }
01715   }
01716   return CUBIT_PNT_BOUNDARY;
01717 }
01718 
01719 
01720 
01721 double CubitOctree::get_scaled_from_wrld_size_grid_node( double wrld_size ){
01722 
01723   double scaled_size;
01724     
01725   scaled_size = 1.0 + ( ( wrld_size - minSizeGridNode ) / ( maxSizeGridNode - minSizeGridNode ) ) * 9.0;
01726 
01727     //PRINT_DEBUG_157("Scaled Size = %2.10f\n",scaled_size);
01728     //PRINT_INFO(" world = %f, scaled = %f\n", wrld_size, scaled_size );
01729   if( scaled_size > 10 ){
01730     PRINT_INFO("Testing ");
01731   }
01732 
01733   return scaled_size;
01734 }
01735 
01736 
01737 
01738 
01739 
01740 CubitStatus CubitOctree::circumcenter(CubitVector &a, CubitVector &b, CubitVector &c, CubitVector &center, double &radius)
01741 {
01742   
01743     // Use coordinates relative to point `a' of the triangle. 
01744   CubitVector vec_ba(a, b);
01745   CubitVector vec_ca(a, c);
01746 
01747     // Squares of lengths of the edges incident to `a'. 
01748   double ba_length = vec_ba.length_squared();
01749   double ca_length = vec_ca.length_squared();
01750 
01751     // Cross product of these edges. 
01752     // (Take your chances with floating-point roundoff.)
01753   CubitVector cross_bc = vec_ba * vec_ca;
01754 
01755     // Calculate the denominator of the formulae. 
01756   double denominator = (cross_bc % cross_bc);
01757   if (denominator == 0.0) {return CUBIT_FAILURE;}
01758   assert(denominator != 0.0); // don't think we need this..
01759 
01760     // Calculate offset (from `a') of circumcenter. 
01761   center  = (ba_length * vec_ca - ca_length * vec_ba) * cross_bc * 0.5;
01762   center /= denominator;
01763 
01764     // radius is length from point `a' to center
01765   radius = center.length();
01766 
01767     // Add point `a' to get global coordinate of center
01768   center += a;
01769 
01770   return CUBIT_SUCCESS;
01771 }
01772 
01773 double CubitOctree::capsule_distance_to_facet(const CubitVector &point, CubitFacet *lp_facet, CubitVector &int_point, CubitBoolean use_projection_only)
01774 {
01775     
01776   CubitVector c = lp_facet->center();
01777   CubitVector n = lp_facet->normal();
01778   double d = (point - c)%n;
01779   CubitVector proj = point - d*n;
01780   CubitVector points[3] = {lp_facet->point(0)->coordinates(), lp_facet->point(1)->coordinates(), lp_facet->point(2)->coordinates()};
01781   
01782   if (use_projection_only || CubitOctreeNode::is_intersection_point_contained_inside_facet(proj, points[0], points[1], points[2]))
01783   {
01784     int_point = proj;
01785     return fabs(d);
01786   }
01787 
01788   double dn[3] = {0,0,0};
01789   CubitVector pn[3];
01790   double t = 0;
01791   int i;
01792 
01793   for (i=0; i < 3; ++i)
01794   {
01795     CubitVector dir = points[(i+1)%3] - points[i];
01796     double length = dir.length();
01797     dir.normalize();
01798     t = ((point - points[i])%dir)/length;
01799     if (t >= 0 && t <= 1)
01800     {
01801       pn[i] = (length*t*dir + points[i]);
01802       dn[i] = (point - (t*length*dir + points[i])).length();
01803     }
01804     else
01805     {
01806       if (t < 0)
01807       {
01808         pn[i] = points[i];
01809         dn[i] = (point - points[i]).length();
01810       }
01811       else if (t > 1)
01812       {
01813         pn[i] = points[(i+1)%3];
01814         dn[i] = (point - points[(i+1)%3]).length();
01815       }
01816         //else {PRINT_INFO("Hmmm.... t is not within bounds: %f!!!!!!!!!!!!!!!!11\n", t);}
01817     }
01818   }
01819   
01820   if (dn[0] <= dn[1] && dn[0] <= dn[2]) {int_point = pn[0]; return dn[0];}
01821   if (dn[1] <= dn[0] && dn[1] <= dn[2]) {int_point = pn[1]; return dn[1];}
01822   if (dn[2] <= dn[1] && dn[2] <= dn[0]) {int_point = pn[2]; return dn[2];}
01823 
01824   else
01825   {
01826       //PRINT_INFO("capsule dist not found!!!!!!!!!!!!!1\n");
01827     return -1;
01828   }
01829 }
01830 
01831 
01832 #ifndef NDEBUG
01833 void CubitOctree::write_sizing_info_file( const char *file_name ){
01834      
01835   int i,j,k;
01836   double size;
01837   
01838   FILE *pof = fopen( file_name, "w" );
01839   if( pof == NULL ){
01840     PRINT_INFO("ERROR: Opening Sizing Output File ");
01841     exit(0);
01842   }
01843      
01844   fprintf( pof, "TFD1\n" );
01845   int num_of_cells = (int)pow(2.0, maxDepth-1 );
01846   fprintf( pof, "%d %d %d\n", num_of_cells, num_of_cells, num_of_cells ); 
01847     //fprintf( pof, "%f %f %f %f %f %f\n", bboxMaxY, bboxMinZ, bboxMinX, bboxMinY, bboxMaxZ, bboxMaxX );
01848   fprintf( pof, "%f %f %f %f %f %f\n", bboxMinX, bboxMinY, bboxMinZ, bboxMaxX, bboxMaxY, bboxMaxZ );
01849 
01850   int l = num_of_cells + 1;
01851   int m = num_of_cells + 1;
01852   int n = num_of_cells + 1;
01853   CubitVector point;
01854 
01855   double grid_size = fabs( bboxMinY - bboxMaxY ) / num_of_cells; 
01856   
01857   for( k = 0; k < n; k++ )
01858   {
01859     for( j = 0; j < m; j++ )
01860     {
01861       for ( i = 0; i < l; i++ )
01862       {
01863         point.x( bboxMinX + grid_size * i );
01864         point.y( bboxMinY + grid_size * j );
01865         point.z( bboxMinZ + grid_size * k );
01866         
01867         size = size_at_a_point( point, OCTREE_SIZE_DEFAULT);    
01868         fprintf( pof, "%f %f %f %f %f %f %f %f %f %f %f %f\n",
01869                  1.0, 0.0, 0.0, size,
01870                  0.0, 1.0, 0.0, size,
01871                  0.0, 0.0, 1.0, size );
01872       }
01873     }
01874   }
01875   
01876   fclose( pof );
01877 }
01878 
01879 void CubitOctree::write_matlab_sizing_info_file( const char *file_name ){
01880      
01881   int i,j,k;
01882   double size;
01883   
01884   CubitOctreeNode *grid_node;
01885 
01886   for (i=0; i < gridNodeVector.size(); ++i)
01887   {
01888     grid_node = gridNodeVector.get_and_step();
01889     if (grid_node->get_color() == CUBIT_GREY_INDEX)
01890     {
01891       grid_node->set_size(0.0, OCTREE_SIZE_DEFAULT);
01892     }
01893   }
01894   
01895   FILE *pof = fopen( file_name, "w" );
01896   if( pof == NULL ){
01897     PRINT_INFO("ERROR: Opening Sizing Output File ");
01898     exit(0);
01899   }
01900      
01901   
01902   int num_of_cells = 50;//(int)pow(2.0, MATLAB_OUTPUT_MAX_DEPTH );
01903     //fprintf( pof, "%f %f %f %f %f %f\n", bboxMinX, bboxMinY, bboxMinZ, bboxMaxX, bboxMaxY, bboxMaxZ );
01904 
01905   int l = num_of_cells + 1;
01906   int m = num_of_cells + 1;
01907   int n = num_of_cells + 1;
01908   CubitVector point;
01909 
01910   double grid_size = fabs( bboxMinY - bboxMaxY ) / num_of_cells; 
01911   
01912   for( k = 0; k < n; k++ ){
01913     for( j = 0; j < m; j++ ){
01914       for ( i = 0; i < l; i++ ){            
01915         point.x( bboxMinX + grid_size * i );
01916         point.y( bboxMinY + grid_size * j );
01917         point.z( bboxMinZ + grid_size * k );
01918         
01919         size = size_at_a_point( point, OCTREE_SIZE_DEFAULT);    
01920 
01921         fprintf( pof, "%f %f %f %f\n",
01922                  point.x(), point.y(), point.z(), size );
01923 
01924       }
01925     }
01926   }
01927   
01928   fclose( pof );
01929 } 
01930 
01931 void CubitOctree::write_octree_sizing_info_file( const char *file_name ){
01932      
01933   int i;
01934   double size;
01935   CubitOctreeNode *ptr_onode;
01936   double FACTOR_FOR_BUBBLEMESH = 1 / 1.8;  // Based on 20% overlap of bubbles
01937   
01938   FILE *pof = fopen( file_name, "w" );
01939   if( pof == NULL ){
01940     PRINT_INFO("ERROR: Opening Sizing Output File ");
01941     exit(0);
01942   }
01943   
01944   DLIList<CubitOctreeCell*> stack;
01945   CubitOctreeCell *ptr_cell;
01946 
01947   fprintf( pof, "OCT 1\n" );
01948 //   int num_of_cells = (int)pow(2.0, skeletonProxy->get_max_depth()-1 );
01949   fprintf( pof, "B %f %f %f %f %f %f\n", 
01950            bboxMinX, bboxMinY, bboxMinZ, bboxMaxX, bboxMaxY, bboxMaxZ );
01951 
01952   gridNodeVector.reset();
01953   for( i = 0; i < gridNodeVector.size(); i++ ){
01954     ptr_onode = gridNodeVector.get_and_step();
01955     size = ptr_onode->get_size(OCTREE_SIZE_DEFAULT);
01956         
01957     fprintf( pof, "V %d %f %f %f %f \n", 
01958              ptr_onode->get_num(), ptr_onode->x(), ptr_onode->y(), 
01959              ptr_onode->z(), size *  FACTOR_FOR_BUBBLEMESH  );
01960   }
01961 
01962   stack.append( root );
01963 
01964   while( stack.size() > 0 ){
01965     ptr_cell = stack.remove();
01966     ptr_cell->write_octreecell_sizing_info_file( pof, stack );
01967   }
01968     
01969   fprintf( pof, "X\n");
01970   fclose( pof );
01971 }
01972 
01973 #endif
01974 
01975 
01976 
01977 //EOF
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines