cgma
ChordalAxis.cpp
Go to the documentation of this file.
00001 //-------------------------------------------------------------------------
00002 // Filename      : ChordalAxis.cpp
00003 //
00004 // Purpose       :
00005 //
00006 // Special Notes : Delaunay Triangulation Using a Uniform Grid
00007 //
00008 // Creator       : Jit Ken Tan
00009 //
00010 // Creation Date : 01/20/03
00011 //-------------------------------------------------------------------------
00012 
00013 #include "ChordalAxis.hpp"
00014 #include "DLIList.hpp"
00015 #include "CubitPoint.hpp"
00016 #include "CubitPointData.hpp"
00017 #include "CubitFacet.hpp"
00018 #include "CubitFacetData.hpp"
00019 #include "CubitFacetEdge.hpp"
00020 #include "CubitFacetEdgeData.hpp"
00021 #include "CubitMessage.hpp"
00022 #include "CpuTimer.hpp"
00023 #include "CubitVector.hpp"
00024 #include "assert.h"
00025 #include "FacetorTool.hpp"
00026 #include "TDChordal.hpp"
00027 #include "ParamTool.hpp"
00028 #include "PlanarParamTool.hpp"
00029 
00030 //-------------------------------------------------------------------------
00031 // Purpose       : Constructor
00032 //
00033 // Special Notes : 
00034 //
00035 // Creator       : Jit Ken Tan
00036 //
00037 // Creation Date : 01/20/03
00038 //-------------------------------------------------------------------------
00039 ChordalAxis::ChordalAxis()
00040 {
00041   startNodes = NULL;
00042   endNodes = NULL;
00043 }
00044 
00045 //-------------------------------------------------------------------------
00046 // Purpose       : Destructor
00047 //
00048 // Special Notes : 
00049 //
00050 // Creator       : Jit Ken Tan
00051 //
00052 // Creation Date : 01/20/03
00053 //-------------------------------------------------------------------------
00054 ChordalAxis::~ChordalAxis()
00055 {
00056   if(startNodes)
00057     delete [] startNodes;
00058   if(endNodes)
00059     delete [] endNodes;
00060 
00061 }
00062 
00063 //-------------------------------------------------------------------------
00064 // Purpose       : Create a delaunay constrained triangulation
00065 //
00066 // Special Notes : 
00067 //
00068 // Creator       : Jit Ken Tan
00069 //
00070 // Creation Date : 01/20/03
00071 //-------------------------------------------------------------------------
00072 
00073 CubitStatus ChordalAxis::create_triangulation(DLIList <CubitPoint*> boundary_loop_list){ 
00074              
00075   int ii;
00076   //int num_boundary_edges = 0;
00077   int dummy_variable = -1;
00078   int surf_id = 1;
00079   CubitStatus status = CUBIT_SUCCESS;
00080   DLIList<CubitPoint*> bounding_nodes;
00081   ParamTool *p_tool = NULL;
00082  
00083   //computing number of boundary edges
00084   numBoundaryEdges = boundary_loop_list.size();
00085 
00086   //finding the start and end nodes of edges
00087   startNodes = new CubitPoint* [numBoundaryEdges];
00088   endNodes = new CubitPoint* [numBoundaryEdges];
00089 
00090   for(ii = boundary_loop_list.size(); ii > 0; ii--){
00091     startNodes[ii-1] = boundary_loop_list.get_and_step();
00092     endNodes[ii-1] = boundary_loop_list.get();
00093   }
00094 
00095    //obtaining a list of boundary_nodes
00096   for(ii = 0; ii < numBoundaryEdges; ii++){
00097     bounding_nodes.append(startNodes[ii]);
00098   }
00099 
00100   FacetorTool<int,CubitFacet, CubitFacetEdge, CubitPoint, CubitFacetData, CubitPointData, int> facetor = FacetorTool<int,CubitFacet, CubitFacetEdge, CubitPoint, CubitFacetData, CubitPointData, int> (&surf_id, bounding_nodes, startNodes, endNodes, numBoundaryEdges, &dummy_variable, p_tool);
00101   
00102   status = facetor.mesh_surfwoIP(triangulation);
00103 
00104   return status;
00105 }
00106 
00107 //-------------------------------------------------------------------------
00108 // Purpose       : Obtain the TDChordal Tool Data
00109 //
00110 // Special Notes : 
00111 //
00112 // Creator       : Jit Ken Tan
00113 //
00114 // Creation Date : 01/20/03
00115 //-------------------------------------------------------------------------
00116 TDChordal * ChordalAxis::get_tool_data(CubitFacet *curr_facet){
00117 
00118   ToolData *td = curr_facet-> get_TD( TDChordal::is_chordal );
00119   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00120   return td_chordal;
00121 }
00122 
00123 //-------------------------------------------------------------------------
00124 // Purpose       : Perform the entire chordal axis extraction
00125 //
00126 // Special Notes : 
00127 //
00128 // Creator       : Jit Ken Tan
00129 //
00130 // Creation Date : 01/20/03
00131 //-------------------------------------------------------------------------
00132 CubitStatus ChordalAxis::chordal_axis_transform(DLIList <CubitPoint*> boundary_loop_list,DLIList <CubitFacet *> &facet_list, DLIList <CubitFacetEdge *> &axis){
00133   CubitStatus status = CUBIT_SUCCESS;
00134 
00135   // initialize alpha
00136   pruningAlpha = 0;
00137   
00138   status = create_triangulation(boundary_loop_list);
00139   if(status != CUBIT_SUCCESS){
00140     PRINT_ERROR("Delaunay Triangulation failed\n");
00141     return status;
00142   }
00143 
00144   facet_list = triangulation;
00145 
00146    status = flagging_boundary();
00147   if(status != CUBIT_SUCCESS){
00148     PRINT_ERROR("Flagging of boundary failed\n");
00149     return status;
00150   }
00151 
00152   status = process_triangles();
00153   if(status != CUBIT_SUCCESS){
00154     PRINT_ERROR("Processing of triangles failed\n");
00155     return status;
00156   }
00157 
00158   if(pruningAlpha > 0){
00159     status = pruning();
00160     if(status != CUBIT_SUCCESS){
00161       PRINT_ERROR("Pruning failed\n");
00162       return status;
00163     }
00164   }
00165 
00166   status =  axis_extraction(axis);
00167   if(status != CUBIT_SUCCESS){
00168     PRINT_ERROR("Axis extraction failed\n");
00169     return status;
00170   }
00171 
00172   axis.reset();
00173 
00174   
00175   return CUBIT_SUCCESS;
00176 }
00177 
00178 //-------------------------------------------------------------------------
00179 // Purpose       : Consolidate the junction triangles in our mesh 
00180 //
00181 // Special Notes : If pruning is implemented, the pruned triangles are ignored //                 in determining junction tris
00182 //
00183 // Creator       : Jit Ken Tan
00184 //
00185 // Creation Date : 01/20/03
00186 //-------------------------------------------------------------------------
00187 CubitStatus ChordalAxis::junction_tris_collection(){
00188   junctionTris.clean_out();
00189   int ii;
00190   for (ii = triangulation.size(); ii > 0; ii--){
00191     CubitFacet *current_facet = triangulation.get_and_step();
00192     TDChordal *td_chordal = get_tool_data(current_facet);
00193     TriType tri_class = td_chordal->get_tritype();
00194 
00195     if(tri_class == JUNCTION)
00196       junctionTris.append(current_facet);
00197   }
00198   return CUBIT_SUCCESS;
00199 }
00200 
00201 //-------------------------------------------------------------------------
00202 // Purpose       : Perform pruning of the chordal axis
00203 //
00204 // Special Notes : 
00205 //
00206 // Creator       : Jit Ken Tan
00207 //
00208 // Creation Date : 01/20/03
00209 //-------------------------------------------------------------------------
00210 CubitStatus ChordalAxis::pruning(){
00211   int ii;
00212   CubitStatus status = CUBIT_SUCCESS;
00213   bool pruning_performed = FALSE;
00214 
00215   junction_tris_collection();
00216   
00217   for(ii = 0; ii < junctionTris.size(); ii++){
00218     CubitFacet *curr_facet = junctionTris.get_and_step();
00219     TDChordal *td_chordal = get_tool_data(curr_facet);
00220     
00221     //error, return failure
00222     if(td_chordal == NULL){
00223       PRINT_ERROR("Tool Data of facet doesn't exist\n");
00224       return CUBIT_FAILURE;
00225     }
00226 
00227     //triangle is already pruned
00228     if(td_chordal->get_tritype() == DISCARDED)
00229       continue;
00230     
00231     int jj;
00232     for(jj = 0; jj < 3; jj++){
00233       mark_visited(curr_facet);
00234       status = junction_triangle_pruning(curr_facet, jj, pruning_performed);
00235       if(status == CUBIT_FAILURE){
00236     PRINT_ERROR("Junction triangle pruning function failed\n");
00237     return status;
00238       }
00239     }
00240 
00241     unmark_visited(curr_facet);
00242 
00243   }
00244 
00245   if(pruning_performed == TRUE){
00246     status = pruning();
00247     if(status == CUBIT_FAILURE)
00248       return status;
00249   }
00250   
00251   return CUBIT_SUCCESS;
00252   
00253 }
00254 
00255 //-------------------------------------------------------------------------
00256 // Purpose       : Pruning commences with a junction triangle
00257 //
00258 // Special Notes : 
00259 //
00260 // Creator       : Jit Ken Tan
00261 //
00262 // Creation Date : 01/20/03
00263 //-------------------------------------------------------------------------
00264 CubitStatus ChordalAxis::junction_triangle_pruning(CubitFacet *curr_facet, int edge_index, bool &pruning_performed){
00265   
00266   CubitStatus status = CUBIT_SUCCESS;
00267   CubitFacetEdge *curr_edge = curr_facet->edge(edge_index);
00268   CubitPoint *start_node = curr_edge->start_node();
00269   CubitPoint *end_node = curr_edge->end_node();
00270   CubitVector start_vector = start_node->coordinates();
00271   CubitVector end_vector = end_node->coordinates();
00272   CubitVector edge_vector = end_vector - start_vector;
00273   double edge_length = edge_vector.length();
00274   CubitVector unit_edge_vector = edge_vector/edge_length;
00275 
00276   CubitFacet *adj_triangle = curr_facet->shared_facet(start_node, end_node);
00277 
00278   if(adj_triangle == NULL){
00279     PRINT_INFO("Error in Collection of Junction Triangles\n");
00280     return CUBIT_FAILURE;
00281   }
00282   
00283   DLIList <CubitFacet *> pruned_triangles;
00284   double pruning_distance = pruningAlpha*edge_length;
00285   CubitStatus prune = CUBIT_SUCCESS;
00286 
00287   status = search_triangles_to_prune(adj_triangle,curr_edge, start_node, unit_edge_vector, pruning_distance, pruned_triangles, prune );
00288   if(status == CUBIT_FAILURE){
00289     PRINT_ERROR("The search for triangles to Prune has failed\n");
00290     return status;
00291   }  
00292 
00293   // We pruned the triangles
00294   if(prune == CUBIT_SUCCESS){
00295     pruning_performed = TRUE;
00296     flag_boundary(curr_facet, start_node,end_node);
00297     determine_tritype(curr_facet);
00298     int kk;
00299     for(kk = pruned_triangles.size(); kk > 0; kk--){
00300       CubitFacet *pruning_triangle = pruned_triangles.get_and_step();
00301       status = mark_tri_discarded(pruning_triangle);
00302       if(status == CUBIT_FAILURE){
00303     PRINT_ERROR("Marking tris as discarded has failed\n");
00304     return CUBIT_FAILURE;
00305       }
00306     }
00307     
00308   }
00309 
00310   // unmarked the visited triangles
00311   int kk;
00312   for(kk = pruned_triangles.size(); kk > 0; kk--){
00313     CubitFacet *visited_triangle = pruned_triangles.get_and_step();
00314     status = unmark_visited(visited_triangle);
00315     if(status == CUBIT_FAILURE){
00316       PRINT_ERROR("Marking tris as unvisited has failed\n");
00317       return CUBIT_FAILURE;
00318     }
00319   }
00320 
00321   return CUBIT_SUCCESS;
00322 
00323 }
00324 
00325 //-------------------------------------------------------------------------
00326 // Purpose       : Marked a triangle as discarded
00327 //
00328 // Special Notes : 
00329 //
00330 // Creator       : Jit Ken Tan
00331 //
00332 // Creation Date : 01/20/03
00333 //-------------------------------------------------------------------------
00334 CubitStatus ChordalAxis::mark_tri_discarded(CubitFacet *curr_facet){
00335 
00336   TDChordal *td_chordal = get_tool_data(curr_facet);
00337   
00338   //error, return failure
00339   if(td_chordal == NULL)
00340     return CUBIT_FAILURE;
00341   
00342   td_chordal->set_tritype(DISCARDED);
00343  
00344   return CUBIT_SUCCESS;
00345 }
00346 
00347 //-------------------------------------------------------------------------
00348 // Purpose       : Searching for triangles to prune
00349 //
00350 // Special Notes : 
00351 //
00352 // Creator       : Jit Ken Tan
00353 //
00354 // Creation Date : 01/20/03
00355 //-------------------------------------------------------------------------
00356 CubitStatus ChordalAxis::search_triangles_to_prune(CubitFacet *curr_facet, CubitFacetEdge *curr_edge, CubitPoint *start_node, CubitVector unit_edge_vector, double pruning_distance, DLIList <CubitFacet *> &pruned_triangles, CubitStatus &prune ){
00357   
00358   CubitStatus status = CUBIT_SUCCESS;
00359   
00360   int ii, other_index;
00361   other_index = curr_facet->other_index(curr_edge->start_node(),curr_edge->end_node());
00362 
00363 
00364   if(other_index == -1){
00365     PRINT_ERROR("Can't find the other_index to the current Facet\n");
00366     return CUBIT_FAILURE;
00367   }
00368 
00369   //Computing perpendicular distance of other point from edge
00370   
00371   CubitVector other_point_vec = curr_facet->point(other_index)->coordinates();
00372   CubitVector edge_to_point_vec = other_point_vec - start_node->coordinates();
00373   double edge_to_point_dist = (edge_to_point_vec * unit_edge_vector).length();
00374   
00375   //make edge_to_point distance positive
00376   if(edge_to_point_dist < 0){
00377     edge_to_point_dist = -edge_to_point_dist;
00378   }
00379  
00380   if(edge_to_point_dist > pruning_distance){
00381     prune = CUBIT_FAILURE;
00382     return CUBIT_SUCCESS;
00383   }
00384 
00385   status = mark_visited(curr_facet);
00386   if(status == CUBIT_FAILURE){
00387     PRINT_ERROR("Marking facet as unvisited has failed\n");
00388     return status;
00389   }
00390 
00391   pruned_triangles.append(curr_facet);
00392 
00393   for(ii = 0; ii < 3 && prune == CUBIT_SUCCESS; ii++){
00394     CubitFacetEdge *curr_edge = curr_facet->edge(ii);
00395     CubitPoint *curr_edge_start = curr_edge->start_node();
00396     CubitPoint *curr_edge_end = curr_edge->end_node();
00397     CubitFacet *adj_facet = curr_facet->shared_facet(curr_edge_start, curr_edge_end);
00398     
00399     // Check if there is an adjacent triangle
00400     if(adj_facet == NULL )
00401       continue;
00402     
00403     // Check if adjacent triangle is visited  
00404     bool visited = TRUE;
00405     status = get_visited(adj_facet, visited);
00406     if(status == CUBIT_FAILURE){
00407       PRINT_ERROR("getting visited status of triangle has failed\n");
00408       return status;
00409     }
00410     if(visited == TRUE)
00411       continue;
00412 
00413     status = search_triangles_to_prune(adj_facet,curr_edge, start_node, unit_edge_vector, pruning_distance, pruned_triangles, prune);
00414     if(status == CUBIT_FAILURE){
00415       PRINT_ERROR("Search for triangles has failed\n");
00416       return CUBIT_FAILURE;
00417     }
00418   }
00419 
00420   return CUBIT_SUCCESS;
00421 }
00422 
00423 //-------------------------------------------------------------------------
00424 // Purpose       : Flagging all the triangles with boundary edges
00425 //
00426 // Special Notes : 
00427 //
00428 // Creator       : Jit Ken Tan
00429 //
00430 // Creation Date : 01/20/03
00431 //-------------------------------------------------------------------------
00432 CubitStatus ChordalAxis::flagging_boundary(){
00433   int ii;
00434   for (ii = 0; ii < numBoundaryEdges; ii++){
00435     CubitFacet *f1, *f2;
00436     startNodes[ii] -> shared_facets(endNodes[ii], f1, f2);
00437     if(!f1 || f2)
00438       return CUBIT_FAILURE;
00439     flag_boundary(f1, startNodes[ii], endNodes[ii]);
00440   }
00441   
00442   return CUBIT_SUCCESS;
00443 
00444 }
00445 
00446 //-------------------------------------------------------------------------
00447 // Purpose       : Mark Triangle as visited
00448 //
00449 // Special Notes : 
00450 //
00451 // Creator       : Jit Ken Tan
00452 //
00453 // Creation Date : 01/20/03
00454 //-------------------------------------------------------------------------
00455 CubitStatus ChordalAxis::mark_visited(CubitFacet *curr_facet){
00456   ToolData *td = curr_facet->get_TD( TDChordal::is_chordal );
00457   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00458 
00459   
00460   if (td_chordal == NULL) {
00461     return CUBIT_FAILURE;
00462   }
00463 
00464   td_chordal->mark_visited();
00465   return CUBIT_SUCCESS;
00466 }
00467 
00468 //-------------------------------------------------------------------------
00469 // Purpose       : Mark Triangle as un-visited
00470 //
00471 // Special Notes : 
00472 //
00473 // Creator       : Jit Ken Tan
00474 //
00475 // Creation Date : 01/20/03
00476 //-------------------------------------------------------------------------
00477 CubitStatus ChordalAxis::unmark_visited(CubitFacet *curr_facet){
00478   ToolData *td = curr_facet->get_TD( TDChordal::is_chordal );
00479   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00480 
00481   
00482   if (td_chordal == NULL) {
00483     return CUBIT_FAILURE;
00484   }
00485   
00486   td_chordal->unmark_visited();
00487   return CUBIT_SUCCESS;
00488 }
00489 
00490 //-------------------------------------------------------------------------
00491 // Purpose       : Check if the triangle has been visited
00492 //
00493 // Special Notes : 
00494 //
00495 // Creator       : Jit Ken Tan
00496 //
00497 // Creation Date : 01/20/03
00498 //-------------------------------------------------------------------------
00499 CubitStatus ChordalAxis::get_visited(CubitFacet *curr_facet, bool &visited){
00500   ToolData *td = curr_facet->get_TD( TDChordal::is_chordal );
00501   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00502 
00503   
00504   if (td_chordal == NULL) {
00505     return CUBIT_FAILURE;
00506   }
00507   
00508   visited = td_chordal->get_visited();
00509   return CUBIT_SUCCESS;
00510 }
00511 
00512 //-------------------------------------------------------------------------
00513 // Purpose       : Mark Triangle with boundary edge
00514 //
00515 // Special Notes : 
00516 //
00517 // Creator       : Jit Ken Tan
00518 //
00519 // Creation Date : 01/20/03
00520 //-------------------------------------------------------------------------
00521 CubitStatus ChordalAxis::flag_boundary(CubitFacet *curr_facet, CubitPoint *start_node, CubitPoint *end_node){
00522   ToolData *td = curr_facet->get_TD( TDChordal::is_chordal );
00523   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00524   int index;
00525 
00526   if (td_chordal == NULL) {
00527     td_chordal = new TDChordal();
00528     curr_facet -> add_TD( td_chordal);
00529   }
00530   
00531   int dummy;
00532   index =  curr_facet->edge_index(start_node, end_node, dummy);
00533   return  (td_chordal->flag_boundary_edge(index));
00534 
00535 }
00536 
00537 //-------------------------------------------------------------------------
00538 // Purpose       : Classifying the triangles based on number of boundary edges
00539 //
00540 // Special Notes : 
00541 //
00542 // Creator       : Jit Ken Tan
00543 //
00544 // Creation Date : 01/20/03
00545 //-------------------------------------------------------------------------
00546 CubitStatus ChordalAxis::process_triangles(){
00547   int ii;
00548   CubitStatus status = CUBIT_SUCCESS;
00549 
00550   for( ii = triangulation.size(); ii > 0; ii--){
00551     CubitFacet *curr_facet = triangulation.get_and_step();
00552     ToolData *td = curr_facet-> get_TD( TDChordal::is_chordal );
00553     TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00554    
00555     // it is a junction triangle if there is no tool attached
00556     if (td_chordal == NULL){
00557       td_chordal = new TDChordal();
00558       curr_facet-> add_TD(td_chordal);
00559     }          
00560     
00561     status = td_chordal -> determine_tritype();
00562     if(status == CUBIT_FAILURE)
00563       return status;
00564  
00565   }
00566   return CUBIT_SUCCESS;
00567 }
00568 
00569 //-------------------------------------------------------------------------
00570 // Purpose       : Determine triangle type based on number of boundary edges
00571 //
00572 // Special Notes : 2 boundary edges -> termination triangle
00573 //                 1 boundary edges -> sleeve triangle
00574 //                 0 boundary edges -> junction triangle
00575 //
00576 // Creator       : Jit Ken Tan
00577 //
00578 // Creation Date : 01/20/03
00579 //-------------------------------------------------------------------------
00580 CubitStatus ChordalAxis::determine_tritype(CubitFacet *curr_facet){
00581   CubitStatus status = CUBIT_SUCCESS;
00582   ToolData *td = curr_facet-> get_TD( TDChordal::is_chordal );
00583   TDChordal *td_chordal = dynamic_cast<TDChordal *> (td);
00584     
00585   if (td_chordal == NULL){
00586     return CUBIT_FAILURE;
00587   }          
00588     
00589   status = td_chordal->determine_tritype();
00590 
00591   return status;
00592 }
00593 
00594 //-------------------------------------------------------------------------
00595 // Purpose       : Perform the chordal axis extraction
00596 //
00597 // Special Notes : 
00598 //
00599 // Creator       : Jit Ken Tan
00600 //
00601 // Creation Date : 01/20/03
00602 //-------------------------------------------------------------------------
00603 CubitStatus ChordalAxis::axis_extraction(DLIList <CubitFacetEdge *> &axis){
00604  
00605   triangulation.reset();
00606   int ii;
00607   TriType tri_class;
00608   CubitStatus status;
00609 
00610   for (ii = triangulation.size(); ii > 0; ii--){
00611     CubitFacet *current_facet = triangulation.get_and_step();
00612     TDChordal *td_chordal = get_tool_data(current_facet);
00613     tri_class = td_chordal->get_tritype();
00614     if(tri_class == UNDEFINED){
00615       PRINT_ERROR("Tri_class is undefined\n");
00616       return CUBIT_FAILURE;
00617     } else if(tri_class == SLEEVE) {
00618 
00619       status = sleeve_triangle(current_facet, td_chordal, axis);
00620       if(status == CUBIT_FAILURE){
00621     PRINT_ERROR("Error in Sleeve Triangles\n");
00622     return status;
00623       }
00624 
00625     } else if(tri_class == JUNCTION){
00626       junction_triangle(current_facet, td_chordal,axis);
00627     }
00628   }
00629   
00630   return CUBIT_SUCCESS;
00631 }
00632 
00633 //-------------------------------------------------------------------------
00634 // Purpose       : Perform chordal axis extraction given a sleeve triangle
00635 //
00636 // Special Notes : 
00637 //
00638 // Creator       : Jit Ken Tan
00639 //
00640 // Creation Date : 01/20/03
00641 //-------------------------------------------------------------------------
00642 CubitStatus ChordalAxis::sleeve_triangle(CubitFacet* curr_facet,TDChordal* td_chordal, DLIList <CubitFacetEdge *> &axis){
00643   
00644   int index1, index2;
00645   DLIList<int> edge_list;
00646   CubitPoint *point11, *point12, *point21, *point22;
00647   CubitVector vector11, vector12, vector21, vector22;
00648   CubitPoint *mid_point1, *mid_point2;
00649   CubitVector mid1, mid2;
00650   CubitFacetEdge *segment;
00651 
00652   //get non_boundary edges
00653   td_chordal->get_non_boundary_edges(edge_list);
00654 
00655   if(edge_list.size() != 2){
00656     PRINT_ERROR("Error in Sleeve Triangles: Boundary Edges Wrong. \n");
00657   }
00658   
00659   index1 = edge_list.get_and_step();
00660   index2 = edge_list.get_and_step();
00661 
00662   curr_facet->get_edge_pts(index1, point11, point12);
00663   curr_facet->get_edge_pts(index2, point21, point22);
00664 
00665   vector11 = point11->coordinates();
00666   vector12 = point12->coordinates();
00667   vector21 = point21->coordinates();
00668   vector22 = point22->coordinates();
00669 
00670   mid1 = (vector11 + vector12)/2;
00671   mid2 = (vector21 + vector22)/2;
00672 
00673   mid_point1 = (CubitPoint *) new CubitPointData(mid1);
00674   mid_point2 = (CubitPoint *) new CubitPointData(mid2);    
00675  
00676   segment = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point1, mid_point2);
00677   
00678   axis.append(segment);
00679   
00680   return CUBIT_SUCCESS;
00681 }
00682 
00683 //-------------------------------------------------------------------------
00684 // Purpose       : Perform chordal axis extraction given a junction triangle
00685 //
00686 // Special Notes : 
00687 //
00688 // Creator       : Jit Ken Tan
00689 //
00690 // Creation Date : 01/20/03
00691 //-------------------------------------------------------------------------
00692 CubitStatus ChordalAxis::junction_triangle(CubitFacet *current_facet, TDChordal* td_chordal, DLIList <CubitFacetEdge *> &axis){
00693  
00694   CubitPoint *point11, *point12, *point21, *point22, *point31, *point32;
00695   CubitVector vector11, vector12, vector21, vector22, vector31, vector32;
00696   double mag_edge_one, mag_edge_two, mag_edge_three;
00697   CubitPoint *mid_point1, *mid_point2, *mid_point3;
00698   CubitVector mid1, mid2, mid3;
00699   CubitVector center;
00700   CubitPoint *center_point;
00701   CubitFacetEdge *segment1 = NULL, *segment2 = NULL, *segment3 = NULL;
00702   CubitVector facet_vec1, facet_vec2, facet_vec3;
00703   bool obtuse = false;
00704   int largest_edge_index = -1;
00705 
00706   // obtain the points of the edges 
00707   current_facet->get_edge_pts(0, point11, point12);
00708   current_facet->get_edge_pts(1, point21, point22);
00709   current_facet->get_edge_pts(2, point31, point32);
00710   
00711   //get vector equivalents of the points
00712   vector11 = point11->coordinates();
00713   vector12 = point12->coordinates();
00714   vector21 = point21->coordinates();
00715   vector22 = point22->coordinates();
00716   vector31 = point31->coordinates();
00717   vector32 = point32->coordinates();
00718 
00719   
00720   //calculate the length of each edge
00721   mag_edge_one = (vector11 - vector12).length();
00722   mag_edge_two = (vector21 - vector22).length();
00723   mag_edge_three = (vector31 - vector32).length();
00724     
00725   //Checking whether triangle is acute or obtuse
00726   if(pow(mag_edge_one,2) >= pow(mag_edge_two,2) + pow(mag_edge_three,2)){
00727     obtuse = true;
00728     largest_edge_index = 1;
00729   } else if(pow(mag_edge_two,2) >= pow(mag_edge_one,2) + pow(mag_edge_three,2)){
00730     obtuse = true;
00731     largest_edge_index = 2;
00732   }else if (pow(mag_edge_three,2) >= pow(mag_edge_one,2) +  pow(mag_edge_two,2)){
00733     obtuse = true;
00734     largest_edge_index = 3;
00735   } else {
00736     obtuse = false;
00737   }
00738 
00739   mid1 =(vector11+vector12)/2;
00740   mid2 =(vector21+vector22)/2;
00741   mid3 =(vector31+vector32)/2;
00742   
00743   mid_point1 = (CubitPoint *) new CubitPointData(mid1);
00744   mid_point2 = (CubitPoint *) new CubitPointData(mid2);
00745   mid_point3 = (CubitPoint *) new CubitPointData(mid3);
00746 
00747   // Extracting the chordal axis.
00748   if(obtuse == false){
00749     //triangle is acute.
00750     //Connect the circumcenter and the midpoints of the edges
00751     facet_vec1 = (current_facet->point(0)) -> coordinates();
00752     facet_vec2 = (current_facet->point(1)) -> coordinates();
00753     facet_vec3 = (current_facet->point(2)) -> coordinates();
00754     
00755     center = get_center(facet_vec1, facet_vec2, facet_vec3);
00756 
00757     center_point = (CubitPoint *) new CubitPointData(center);
00758 
00759     segment1 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point1, center_point);
00760     segment2 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point2, center_point);
00761     segment3 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point3, center_point);
00762 
00763     axis.append(segment1);
00764     axis.append(segment2);
00765     axis.append(segment3);
00766 
00767   
00768   } else {
00769     //triangle is obtuse
00770 
00771     switch(largest_edge_index){
00772     case 1:
00773       segment1 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point1, mid_point2);
00774       segment2 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point1, mid_point3);
00775       break;
00776     case 2:
00777       segment1 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point2, mid_point1);
00778       segment2 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point2, mid_point3);
00779       break;
00780     case 3:
00781       segment1 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point3, mid_point1);
00782       segment2 = (CubitFacetEdge *) new CubitFacetEdgeData(mid_point3, mid_point2);
00783       break;
00784     default:
00785       //Not a possible scenario
00786       assert(0); 
00787       break;
00788     }
00789     
00790     axis.append(segment1);
00791     axis.append(segment2);
00792   }
00793 
00794   return CUBIT_SUCCESS;
00795 
00796 }
00797 
00798 //-------------------------------------------------------------------------
00799 // Purpose       : get circumcenter given 3 points of a circle
00800 //
00801 // Special Notes : 
00802 //
00803 // Creator       : Jit Ken Tan
00804 //
00805 // Creation Date : 01/20/03
00806 //-------------------------------------------------------------------------
00807 CubitVector ChordalAxis::get_center(CubitVector v1, CubitVector v2, CubitVector v3){
00808   double x1 = v1.x();
00809   double y1 = v1.y();
00810   double x2 = v2.x();
00811   double y2 = v2.y();
00812   double x3 = v3.x();
00813   double y3 = v3.y();
00814 
00815   double ma = 0.0;
00816   double mb = 0.0;
00817   double cx = 0.0;
00818   double cy = 0.0;
00819 
00820 #ifndef NDEBUG
00821   //for debugging
00822   bool center_found = false;
00823 #endif
00824 
00825   if(x2 != x1 && x3 != x2 ){
00826     ma = (y2-y1)/(x2-x1);
00827     mb = (y3-y2)/(x3-x2);
00828     cx = (ma*mb*(y1-y3) + mb*(x1+x2) - ma*(x2+x3))/(2*(mb-ma));
00829     
00830     if(ma != 0)
00831       cy = (-1/ma)*(cx - (x1+x2)/2)+ (y1+y2)/2;
00832     else
00833       cy = (-1/mb)*(cx - (x3+x2)/2)+ (y3+y2)/2;
00834 
00835 #ifndef NDEBUG
00836     center_found = true;
00837 #endif
00838   } else if (x3 != x1 && x2 != x3) {
00839     ma = (y3-y1)/(x3-x1);
00840     mb = (y2-y3)/(x2-x3);
00841     cx = (ma*mb*(y1-y2) + mb*(x1+x3) - ma*(x2+x3))/(2*(mb-ma));
00842     if(ma != 0)
00843       cy = (-1/ma)*(cx - (x1+x3)/2)+ (y1+y3)/2;
00844     else 
00845       cy = (-1/mb)*(cx - (x3+x2)/2)+ (y3+y2)/2;
00846 
00847 #ifndef NDEBUG
00848     center_found = true;
00849 #endif
00850   } else if ( x1 != x2 && x1 != x3 ){
00851     ma = (y1-y2)/(x1-x2);
00852     mb = (y3-y1)/(x3-x1);    
00853     cx = (ma*mb*(y2-y3) + mb*(x1+x2) - ma*(x1+x3))/(2*(mb-ma));
00854     if(ma != 0)
00855       cy = (-1/ma)*(cx - (x1+x2)/2)+ (y1+y2)/2;
00856     else
00857       cy = (-1/mb)*(cx - (x1+x3)/2)+ (y1+y3)/2;
00858 
00859 #ifndef NDEBUG
00860     center_found = true;
00861 #endif
00862   } else {
00863     //scenario not possible
00864     assert(0);
00865   }
00866 
00867   assert(center_found);
00868 
00869   return CubitVector(cx,cy,0);
00870 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines