cgma
FacetorUtil.cpp
Go to the documentation of this file.
00001 
00002 #ifdef INLINE_TEMPLATES
00003 #define MY_INLINE inline
00004 #else
00005 #define MY_INLINE
00006 #endif
00007 
00008 #include "FacetorUtil.hpp"
00009 #include "TDDelaunay.hpp"
00010 #include "FacetEvalTool.hpp"
00011 #include "CubitPoint.hpp"
00012 #include "CubitPointData.hpp"
00013 #include "CubitFacet.hpp"
00014 #include "CubitFacetData.hpp"
00015 #include "ParamTool.hpp"
00016 #include "CubitMessage.hpp"
00017 
00018 #define DETERM(p1,q1,p2,q2,p3,q3)\
00019      ((q3)*((p2)-(p1)) + (q2)*((p1)-(p3)) + (q1)*((p3)-(p2)))
00020 #define FT_INSIDE_TOL 1.0e-6
00021 #define SQR(x) ((x) * (x))
00022 #define ALPHA 0.70228615
00023 
00024 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00025 CubitStatus
00026 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::locate_point(
00027                                       CubitVector& the_point,
00028                                       DLIList<TRI*>& facet_list,
00029                                       TRI* starting_tri,
00030                                       SURF* owning_surface,
00031                                       TRI*& tri_ptr)
00032 {
00033   CubitStatus rv = CUBIT_SUCCESS;
00034 
00035     // start with the last one found
00036 
00037   if (starting_tri != NULL)
00038       tri_ptr = starting_tri;
00039 
00040     // otherwise use the first one on the list
00041 
00042   else
00043   {
00044     tri_ptr = facet_list.get();
00045   }
00046 
00047 
00048     // loop until we find something
00049 
00050   NODE *n0, *n1, *n2;
00051   double aa, bb, cc;
00052   CubitBoolean found = CUBIT_FALSE;
00053 
00054     //avoiding infinite loop
00055   int counter = 0;
00056   int max_count = facet_list.size();
00057   
00058   while(!found && rv == CUBIT_SUCCESS)
00059   {
00060     tri_ptr->tri_nodes( n0, n1, n2 );
00061     aa = DETERM(the_point.x(), the_point.y(),
00062                 n1->coordinates().x(), n1->coordinates().y(),
00063                 n2->coordinates().x(), n2->coordinates().y());
00064     bb = DETERM(n0->coordinates().x(), n0->coordinates().y(),
00065                 the_point.x(), the_point.y(),
00066                 n2->coordinates().x(), n2->coordinates().y());
00067     cc = DETERM(n0->coordinates().x(), n0->coordinates().y(),
00068                 n1->coordinates().x(), n1->coordinates().y(),
00069                 the_point.x(), the_point.y());
00070     if (aa > -FT_INSIDE_TOL &&
00071         bb > -FT_INSIDE_TOL &&
00072         cc > -FT_INSIDE_TOL)
00073     {
00074       found = CUBIT_TRUE;  // this is the one
00075     }
00076     else
00077     {
00078         // set up to check the next logical neighbor
00079       if (aa <= bb && aa <= cc) 
00080       {
00081         int edge_index = 1;
00082         tri_ptr = tri_ptr->adjacent( edge_index, owning_surface );
00083       }
00084       else if (bb <= aa && bb <= cc) 
00085       {
00086         int edge_index = 2;
00087         tri_ptr = tri_ptr->adjacent( edge_index, owning_surface );
00088       }
00089       else 
00090       {
00091         int edge_index = 0;
00092         tri_ptr = tri_ptr->adjacent( edge_index, owning_surface );
00093       }
00094         // check to see if we've left the triangulation
00095         // also make sure that we are not stuck in a cycle
00096       if (tri_ptr == NULL || counter > max_count)
00097       {
00098         if(counter>max_count){
00099           PRINT_WARNING("Encountered problem locating a triangle; going to exhaustive search.\n");
00100         }
00101 
00102         rv = exhaustive_locate_point( the_point, facet_list, tri_ptr );
00103         found = CUBIT_TRUE;
00104       }
00105     }
00106     ++counter;
00107   }
00108   
00109   return rv;
00110 }
00111 
00112 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00113 CubitStatus
00114 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::exhaustive_locate_point(
00115                                                  CubitVector& the_point,
00116                                                  DLIList<TRI*>& facet_list,
00117                                                  TRI*& tri_ptr)
00118 {
00119   CubitStatus rv = CUBIT_SUCCESS;
00120   
00121   // loop until we find something
00122 
00123   int ii;
00124   NODE *n0, *n1, *n2;
00125   double aa, bb, cc;
00126   CubitBoolean found = CUBIT_FALSE;
00127   for (ii=0; ii<facet_list.size() && !found; ii++)
00128   {
00129     tri_ptr = facet_list.get_and_step();
00130     tri_ptr->tri_nodes( n0, n1, n2 );
00131     aa = DETERM(the_point.x(), the_point.y(),
00132                 n1->coordinates().x(), n1->coordinates().y(),
00133                 n2->coordinates().x(), n2->coordinates().y());
00134     bb = DETERM(n0->coordinates().x(), n0->coordinates().y(),
00135                 the_point.x(), the_point.y(),
00136                 n2->coordinates().x(), n2->coordinates().y());
00137     cc = DETERM(n0->coordinates().x(), n0->coordinates().y(),
00138                 n1->coordinates().x(), n1->coordinates().y(),
00139                 the_point.x(), the_point.y());
00140     if (aa > -FT_INSIDE_TOL &&
00141         bb > -FT_INSIDE_TOL &&
00142         cc > -FT_INSIDE_TOL)
00143     {
00144       found = CUBIT_TRUE;  // this is the one
00145     }
00146   }
00147   if (!found)
00148   {
00149     rv = CUBIT_FAILURE;
00150     tri_ptr = NULL;
00151   }
00152 
00153   return rv;
00154 }
00155 
00156 //-------------------------------------------------------------------------
00157 // Function:    are_nodes_colinear
00158 // Description: determine if the TRI is valid
00159 // Author:      mbrewer
00160 // Date:        6/3/2002
00161 //-------------------------------------------------------------------------
00162 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00163 CubitBoolean
00164 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::are_nodes_colinear(TRI* tri_ptr)
00165 {
00166   NODE *nodes[3];
00167   tri_ptr->tri_nodes( nodes[0], nodes[1], nodes[2] );
00168   double det = DETERM( nodes[0]->coordinates().x(),
00169                        nodes[0]->coordinates().y(),
00170                        nodes[1]->coordinates().x(),
00171                        nodes[1]->coordinates().y(),
00172                        nodes[2]->coordinates().x(),
00173                        nodes[2]->coordinates().y());
00174     //PRINT_INFO("Det = %f\n",det);
00175   
00176   if(fabs(det) > CUBIT_DBL_MIN){
00177     return CUBIT_TRUE;
00178   }
00179   return CUBIT_FALSE;
00180 }
00181 
00182 //-------------------------------------------------------------------------
00183 // Function:    circumcenter
00184 // Description: get the circumcenter of the triangle
00185 // Author:      chynes
00186 // Date:        6/6/2002
00187 //-------------------------------------------------------------------------
00188 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00189 CubitVector& 
00190 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::circumcenter(TRI* tri_ptr)
00191 {
00192     ToolData *td = tri_ptr->get_TD( TDDelaunay<TRI,NODE>::is_delaunay );
00193     TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00194     if(td_del == NULL) {
00195         td_del = new TDDelaunay<TRI, NODE>();
00196         tri_ptr->add_TD( td_del );
00197     }
00198     return td_del->circumcenter2d( tri_ptr );
00199 }
00200 
00201 //-------------------------------------------------------------------------
00202 // Function:    tri_visited
00203 // Description: set the tri_visited flag
00204 // Author:      chynes
00205 // Date:        6/6/2002
00206 //-------------------------------------------------------------------------
00207 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00208 CubitBoolean
00209 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::tri_visited(
00210                                       TRI *tri_ptr,
00211                                       int curr_visit_flag)
00212 {
00213   ToolData *td = tri_ptr->get_TD( TDDelaunay< TRI, NODE >::is_delaunay );
00214   TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00215   if (td_del == NULL)
00216   {
00217     td_del = new TDDelaunay< TRI, NODE >();
00218     tri_ptr->add_TD( td_del );
00219   }
00220   return (td_del->visit_flag() == curr_visit_flag);
00221 }
00222 
00223 //-------------------------------------------------------------------------
00224 // Function:    tri_visited
00225 // Description: set the tri_visited flag
00226 // Author:      chynes
00227 // Date:        6/3/2002
00228 //-------------------------------------------------------------------------
00229 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00230 void
00231 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::tri_visited(
00232                                TRI *facet_ptr,
00233                                CubitBoolean visited,
00234                                int curr_visit_flag)
00235 {
00236   ToolData *td = facet_ptr->get_TD( TDDelaunay< TRI, NODE >::is_delaunay );
00237   TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00238   if (td_del == NULL)
00239   {
00240     td_del = new TDDelaunay< TRI, NODE >();
00241     facet_ptr->add_TD( td_del );
00242   }
00243   if (visited)
00244     td_del->visit_flag(curr_visit_flag);
00245   else
00246     td_del->visit_flag(INT_MIN);
00247 }
00248 
00249 //-------------------------------------------------------------------------
00250 // Function:    radius
00251 // Description: get the radius squared of the triangle circumcircle
00252 // Author:      chynes
00253 // Date:        6/6/2002
00254 //-------------------------------------------------------------------------
00255 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00256 double
00257 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::radius(TRI* tri_ptr)
00258 {
00259   ToolData *td = tri_ptr->get_TD( TDDelaunay< TRI, NODE >::is_delaunay );
00260   TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00261   if (td_del == NULL)
00262   {
00263     td_del = new TDDelaunay< TRI, NODE >();
00264     tri_ptr->add_TD( td_del );
00265   }
00266   return td_del->radius2d( tri_ptr );
00267 }
00268 
00269 //-------------------------------------------------------------------------
00270 // Function:    point_in_circumcircle
00271 // Description: determine if the point is inside the circumcircle of the
00272 //              triangle and recurse to the adjacent triangles
00273 // Author:      chynes
00274 // Date:        6/3/2002
00275 //-------------------------------------------------------------------------
00276 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00277 CubitStatus
00278 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::point_in_circumcircle(
00279                                                CubitVector& the_point,
00280                                                TRI* tri_ptr,
00281                                                DLIList<TRI*>& tri_list,
00282                                                SURF* surface_ptr,
00283                                                int curr_visit_flag)
00284 {
00285   CubitStatus rv = CUBIT_SUCCESS;
00286 
00287   // check this triangle.  If the nodes are colinear do not try to calculate
00288     //the circumcenter.
00289   if(!are_nodes_colinear(tri_ptr))
00290   {
00291     PRINT_ERROR("Can't evaluate center of circumcircle\n");
00292     return CUBIT_FAILURE;
00293   }
00294      
00295   CubitVector cc = circumcenter( tri_ptr );
00296   tri_visited( tri_ptr, CUBIT_TRUE, curr_visit_flag );
00297   double dist2 = SQR(the_point.x() - cc.x()) + SQR(the_point.y() - cc.y());
00298   double r2 = radius( tri_ptr );
00299   double tol_factor = CUBIT_MAX(CUBIT_MAX(tri_ptr->edge(0)->length(),
00300                                           tri_ptr->edge(1)->length()),
00301                                 tri_ptr->edge(2)->length());
00302     //PRINT_INFO("Tolerance factor = %f\n", tol_factor);
00303   
00304   
00305   
00306   if (r2-dist2 > -(tol_factor*FT_INSIDE_TOL*FT_INSIDE_TOL))// inside or on circle
00307   {
00308     tri_list.append( tri_ptr );
00309 
00310     // go to its neighbors
00311 
00312     int iedge;
00313     TRI *adjtri_ptr;
00314     for (iedge=0; iedge<3 && rv == CUBIT_SUCCESS; iedge++){
00315      
00316       int ii = iedge;
00317       adjtri_ptr = tri_ptr->adjacent( ii, surface_ptr );
00318       if (adjtri_ptr != NULL && !tri_visited( adjtri_ptr, curr_visit_flag ))
00319       {
00320         rv = point_in_circumcircle( the_point, adjtri_ptr, tri_list,
00321                                                    surface_ptr, curr_visit_flag );
00322       }
00323     }
00324   }
00325   return rv;
00326 }
00327 
00328 //-------------------------------------------------------------------------
00329 // Function:    natural_neighbor_tris
00330 // Description: get a list of all triangles whose circumcircle contain 
00331 //              the point
00332 // Author:      chynes
00333 // Date:        6/3/2002
00334 //-------------------------------------------------------------------------
00335 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00336 CubitStatus
00337 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::natural_neighbor_tris(
00338                                                CubitVector& the_point,
00339                                                DLIList<TRI*>& facet_list,
00340                                                TRI*& start_tri,
00341                                                SURF* surface_ptr,
00342                                                int& curr_visit_flag,
00343                                                DLIList<TRI*>& tri_list)
00344 {
00345   CubitStatus rv = CUBIT_SUCCESS;
00346 
00347   // find the triangle the point is contained in
00348 
00349   //CubitVector areacoord;
00350   TRI *tri_ptr;
00351   rv = locate_point( the_point, facet_list, start_tri, surface_ptr, tri_ptr );
00352   start_tri = tri_ptr;
00353   
00354   // keep track of visitation to triangles by incrementing curr_visit_flag
00355   // and comparing with the visit flag stored with the triangle
00356 
00357   curr_visit_flag++;
00358 
00359   // Recursively search, (starting with the tri the point is in)
00360   // search for all tris whose circumcircle contain the point and place 
00361   // in the tri_list 
00362 
00363   if (rv == CUBIT_SUCCESS)
00364   {
00365     tri_list.append( tri_ptr );
00366     tri_visited( tri_ptr, CUBIT_TRUE, curr_visit_flag );
00367     int iedge;
00368     TRI *adjtri_ptr;
00369     for (iedge=0; iedge<3 && rv == CUBIT_SUCCESS; iedge++)
00370     {
00371       int ii = iedge;
00372       adjtri_ptr = tri_ptr->adjacent( ii, surface_ptr );
00373       if (adjtri_ptr != NULL && !tri_visited( adjtri_ptr, curr_visit_flag )){
00374         rv = point_in_circumcircle( the_point, adjtri_ptr,
00375                                                    tri_list, surface_ptr,
00376                                                    curr_visit_flag);
00377       }
00378     }
00379   }
00380   return rv;
00381 }
00382 
00383 /************************************************************/
00384 //author:  mbrewer
00385 //This function performs tests on the void created by a point's
00386 //insertion to make sure that invalid connectivities are not
00387 //being created.  
00388 /************************************************************/
00389 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00390 CubitStatus
00391 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::valid_void(
00392                                      NODE * /*point_ptr*/,
00393                                      DLIList<TRI *> &tri_list,
00394                                      SURF* surface_ptr,
00395                                      int curr_visit_flag)
00396 {
00397   int temp_i, temp_j;
00398   DLIList<EDGE*> boundary_edge_list;
00399   DLIList<NODE*> boundary_node_list;
00400   TRI *adjtri_ptr;
00401   TRI *tri_ptr;
00402   EDGE *edge_ptr;
00403   DLIList<EDGE *> edge_list;
00404     //loop over the tri's in tri_list and find all of the curves
00405     //on the boundary of the set (ie, on the boundary of the void). 
00406   for (temp_i=0; temp_i<tri_list.size(); temp_i++)
00407   {
00408     tri_ptr = tri_list.get_and_step();
00409       //check each edge to see whether it is a boundary edge or not
00410     for (temp_j=0; temp_j<3; temp_j++){
00411       
00412       int kk = temp_j;
00413         // - if TRI == CubitFacet
00414         // - kk will be corrected in adjacent() to 
00415         // - give the correct EDGE index
00416       adjtri_ptr = tri_ptr->adjacent( kk, surface_ptr );
00417       if (!adjtri_ptr || !tri_visited( adjtri_ptr, curr_visit_flag ))
00418       {
00419         edge_ptr = tri_ptr->edge(kk);
00420         boundary_edge_list.append(edge_ptr);
00421       }
00422     }
00423   }
00424   int list_size = boundary_edge_list.size();
00425     //uniquify the boundary edge list
00426   boundary_edge_list.uniquify_ordered();
00427     //the list shouldn't have changed size during the uniquify.
00428     //if it did, there is a problem.
00429   if(list_size != boundary_edge_list.size()){
00430     PRINT_WARNING("Unexpected result.  Edge was duplicated on boundary.\n");
00431     return CUBIT_FAILURE;
00432   }
00433     //loop over the boundary edges and get the nodes in the boundary loop
00434   for(temp_i=0; temp_i<list_size; ++temp_i){
00435     edge_ptr=boundary_edge_list.get_and_step();
00436     boundary_node_list.append(edge_ptr->end_node());
00437     boundary_node_list.append(edge_ptr->start_node());
00438   }
00439   list_size = boundary_node_list.size();
00440     //each node should be in exactly two edges.  First make sure that there
00441     //are the correct number of nodes.
00442   if(list_size%2){
00443     PRINT_WARNING("Unexpected result.  Node not listed twice.\n");
00444     return CUBIT_FAILURE;
00445   }
00446     //now uniquify to make sure that the nodes were listed the correct number
00447     //of times.  
00448   boundary_node_list.uniquify_ordered();
00449   if( (list_size/2) != boundary_node_list.size()){
00450     PRINT_WARNING("Unexpected result.  Node was listed an incorrect number of times.\n");
00451     return CUBIT_FAILURE;
00452   }
00453   return CUBIT_SUCCESS;
00454   
00455 }
00456 
00457 //-------------------------------------------------------------------------
00458 // Function:    bowyer_watson_insert
00459 // Description: Bowyer-Watson insertion into an existing Delaunay Mesh
00460 // Author:      chynes
00461 // Date:        6/3/2002
00462 //-------------------------------------------------------------------------
00463 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00464 CubitStatus
00465 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::bowyer_watson_insert(
00466                                               NODE* point_ptr,
00467                                               DLIList<TRI*>& tri_list,
00468                                               DLIList<TRI*>& facet_list,
00469                                               int& curr_visit_flag,
00470                                               SURF* surface_ptr,
00471                                               TRI*& last_tri)
00472 {
00473   CubitStatus rv = CUBIT_SUCCESS;
00474 
00475     // mark the tris in the list so we can distinguish them from their 
00476     // neighbors
00477 
00478   curr_visit_flag++;
00479   int ii, jj;
00480   TRI *tri_ptr;
00481   for (ii=0; ii<tri_list.size(); ii++)
00482   {
00483     tri_ptr = tri_list.get_and_step();
00484     tri_visited( tri_ptr, CUBIT_TRUE, curr_visit_flag );
00485   }
00486     //MBREWER:  This is not an optimal test.  But, we need need to 
00487     //do some tests to try to ensure that the void is valid for what
00488     //we need.  This is attempting to avoid crashes by not allowing nodes
00489     //to be inserted when the mesh starts diverging from the Delaunay
00490     //criteria.
00491   rv = valid_void( point_ptr, tri_list, surface_ptr, curr_visit_flag );
00492   if(!rv)
00493       return rv;
00494   
00495     // find all edges at the boundary of the visited triangles and create
00496     // new triangles with them
00497 
00498     // create a new triangle with this edge and the node
00499   TRI *adjtri_ptr;
00500   TRI *new_tri = NULL;
00501   EDGE *edge_ptr;
00502   DLIList<EDGE *> edge_list;
00503   for (ii=0; ii<tri_list.size(); ii++)
00504   {
00505     tri_ptr = tri_list.get_and_step();
00506     for (jj=0; jj<3; jj++){
00507       
00508       int kk = jj;
00509         // - if TRI == CubitFacet
00510         // - kk will be corrected in adjacent() to 
00511         // - give the correct EDGE index
00512       adjtri_ptr = tri_ptr->adjacent( kk, surface_ptr );
00513       if (!adjtri_ptr || !tri_visited( adjtri_ptr, curr_visit_flag ))
00514       {
00515         edge_ptr = tri_ptr->edge(kk);
00516         assert(edge_list.append_unique( edge_ptr ));
00517         if(tri_ptr->sense(kk) == CUBIT_FORWARD)
00518             new_tri = (TRI *) new TRICHILD( edge_ptr->start_node(), edge_ptr->end_node(), point_ptr, surface_ptr);
00519         else
00520             new_tri = (TRI *) new TRICHILD( edge_ptr->end_node(), edge_ptr->start_node(), point_ptr, surface_ptr);
00521         facet_list.append(new_tri);
00522       }
00523     }
00524   }
00525   last_tri = new_tri;
00526 
00527     // delete the triangles in the original triangle list
00528 
00529   EDGE *del_edge_ptr;
00530   DLIList<EDGE *> del_edge_list;
00531   for (ii=0; ii<tri_list.size(); ii++)
00532   {
00533     tri_ptr = tri_list.get_and_step();
00534     del_edge_list.clean_out();
00535     facet_list.move_to(tri_ptr);
00536     facet_list.extract();
00537     tri_ptr->edges( del_edge_list );
00538     delete tri_ptr;
00539 
00540       // delete the unused edges
00541     for (jj=0; jj<del_edge_list.size(); jj++)
00542     {
00543       del_edge_ptr = del_edge_list.get_and_step();
00544       if (del_edge_ptr->number_tris() == 0 && del_edge_ptr->number_faces() == 0 )
00545           delete del_edge_ptr;
00546     }
00547   }
00548 
00549   return rv;
00550 }
00551 
00552 //-------------------------------------------------------------------------
00553 // Function:    insert_node
00554 // Description: insert one node into Delaunay mesh
00555 // Author:      chynes
00556 // Date:        6/3/2002
00557 //-------------------------------------------------------------------------
00558 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00559 CubitStatus
00560 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::insert_node(
00561                                      NODE *node_ptr,
00562                                      DLIList<TRI*>& facet_list,
00563                                      SURF* surface_ptr,
00564                                      int& curr_visit_flag,
00565                                      TRI*& last_tri)
00566 {
00567   CubitStatus rv = CUBIT_SUCCESS;
00568 
00569   // get a list of all triangles whose circumcircle contain the point
00570   
00571   DLIList<TRI *> tri_list;
00572   CubitVector the_point = node_ptr->coordinates();
00573   rv = natural_neighbor_tris( the_point, facet_list,
00574                               last_tri, surface_ptr,
00575                               curr_visit_flag, tri_list );
00576 
00577 
00578   // Use a Bowyer-Watson insertion 
00579   
00580   if (rv == CUBIT_SUCCESS)
00581   {
00582     rv = bowyer_watson_insert( node_ptr, tri_list,
00583                                facet_list, curr_visit_flag,
00584                                surface_ptr, last_tri);
00585   }
00586 
00587   return rv;
00588 }
00589 
00590 //-------------------------------------------------------------------------
00591 // Function:    get_size
00592 // Description: get the distortion factor for point inside tri, if one exists
00593 //              otherwise return 1;
00594 // Author:      chynes
00595 // Date:        7/24/02
00596 //-------------------------------------------------------------------------
00597 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00598 double
00599 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::get_size(CubitVector &cc, TRI *tri_ptr)
00600 {
00601     //extract data
00602     NODE *n0,*n1,*n2;
00603     CubitVector area_coord;
00604     tri_ptr->tri_nodes(n0,n1,n2);
00605     
00606     if    (n0->coordinates().z() - 1.0 < fabs(FT_INSIDE_TOL) 
00607         && n1->coordinates().z() - 1.0 < fabs(FT_INSIDE_TOL)
00608         && n2->coordinates().z() - 1.0 < fabs(FT_INSIDE_TOL)  )
00609         return 1.0;
00610     else 
00611     {
00612         //get vectors
00613         CubitVector v0 = n0->coordinates();
00614         CubitVector v1 = n1->coordinates();
00615         CubitVector v2 = n2->coordinates();
00616 
00617         //set z direction
00618         v0.z(cc.z());
00619         v1.z(cc.z());
00620         v2.z(cc.z());
00621 
00622         //create points
00623         CubitPoint *p0 = (CubitPoint*) new CubitPointData(v0);
00624         CubitPoint *p1 = (CubitPoint*) new CubitPointData(v1);
00625         CubitPoint *p2 = (CubitPoint*) new CubitPointData(v2);
00626 
00627         //create facet
00628         CubitFacet *temp_facet = (CubitFacet*) new CubitFacetData(p0,p1,p2);
00629 
00630         FacetEvalTool::facet_area_coordinate(temp_facet, cc, area_coord);
00631         
00632         delete p0;
00633         delete p1;
00634         delete p2;
00635         delete temp_facet;
00636 
00637         return    (area_coord.x()*n0->coordinates().z()) 
00638                 + (area_coord.y()*n1->coordinates().z()) 
00639                 + (area_coord.z()*n2->coordinates().z());
00640     }
00641 }
00642 
00643 //-------------------------------------------------------------------------
00644 // Function:    tri_sort_list
00645 // Description: set the tri sort list index
00646 // Author:      chynes
00647 // Date:        6/3/2002
00648 //-------------------------------------------------------------------------
00649 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00650 void
00651 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::tri_sort_list(
00652                                 TRI *facet_ptr,
00653                                 int sort_list_index )
00654 {
00655   ToolData *td = facet_ptr->get_TD( TDDelaunay< TRI, NODE >::is_delaunay );
00656   TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00657   if (td_del == NULL)
00658   {
00659     td_del = new TDDelaunay<TRI, NODE>();
00660     facet_ptr->add_TD( td_del );
00661   }
00662   td_del->tri_sort_list(sort_list_index);
00663 }
00664 
00665 
00666 //-------------------------------------------------------------------------
00667 // Function:    tri_sort_list
00668 // Description: get the tri sort list index
00669 // Author:      chynes
00670 // Date:        6/3/2002
00671 //-------------------------------------------------------------------------
00672 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00673 int
00674 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::tri_sort_list( TRI *facet_ptr )
00675 {
00676   ToolData *td = facet_ptr->get_TD( TDDelaunay< TRI, NODE >::is_delaunay );
00677   TDDelaunay< TRI, NODE > *td_del = dynamic_cast<TDDelaunay< TRI, NODE >*> (td);
00678   if (td_del == NULL)
00679   {
00680     td_del = new TDDelaunay<TRI, NODE>();
00681     facet_ptr->add_TD( td_del );
00682   }
00683   return td_del->tri_sort_list();
00684 }
00685 
00686 //-------------------------------------------------------------------------
00687 // Function:    classify_tri_by_angle
00688 // Description: compute the angles at the triangle vertices and classify
00689 //              by its worst triangle
00690 // Author:      chynes
00691 // Date:        6/3/2002
00692 //-------------------------------------------------------------------------
00693 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00694 CubitStatus
00695 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::classify_tri_by_angle(
00696                                                TRI* tri_ptr,
00697                                                DLIList<TRI*>* sorted_lists,
00698                                                const int num_lists,
00699                                                const double interval,
00700                                                const double quality_angle)
00701 {
00702   //CubitStatus rv = CUBIT_SUCCESS;
00703 
00704   // Determine the minimum angle
00705 
00706   NODE *nodes[3];
00707   tri_ptr->tri_nodes( nodes[0], nodes[1], nodes[2] );
00708   double x0 = nodes[0]->coordinates().x();
00709   double y0 = nodes[0]->coordinates().y();
00710   double x1 = nodes[1]->coordinates().x();
00711   double y1 = nodes[1]->coordinates().y();
00712   double x2 = nodes[2]->coordinates().x();
00713   double y2 = nodes[2]->coordinates().y();
00714 
00715   double ax = x1 - x0;
00716   double ay = y1 - y0;
00717   double bx = x2 - x0;
00718   double by = y2 - y0;
00719   double dot = ax*bx + ay*by;
00720   double a_mag = sqrt( ax*ax + ay*ay );
00721   double b_mag = sqrt( bx*bx + by*by );
00722   double angle0 = dot / ( a_mag * b_mag );
00723   angle0 = acos( angle0 );
00724 
00725   ax = -ax;
00726   ay = -ay;
00727   bx = x2 - x1;
00728   by = y2 - y1;
00729   dot = ax*bx + ay*by;
00730   b_mag = sqrt( bx*bx + by*by );
00731   double angle1 = dot / ( a_mag * b_mag );
00732   angle1 = acos( angle1 );
00733 
00734   double angle2 = CUBIT_PI - angle0 - angle1;
00735 
00736   double min_angle = CUBIT_MIN( CUBIT_MIN(angle0,angle1), 
00737                                 CUBIT_MIN(angle1,angle2) );
00738   if (min_angle < 0.0) {
00739     assert(0);
00740     return CUBIT_FAILURE;
00741   }
00742 
00743   // If it is greater than the QUALITY_ANGLE then place it in
00744   // the triSortArray[0], otherwise place it in one of the other deques
00745   // depending upon its minimum angle
00746 
00747   // Determine which list
00748 
00749   int index;
00750   if (min_angle >= quality_angle) {
00751     index = 0;
00752   }
00753   else {
00754     index = 1 + (int)(min_angle/interval);
00755     if (index < 1) index = 1;
00756     if (index > num_lists-1) index = num_lists-1;
00757   }
00758 
00759   // Place it on a list
00760 
00761   sorted_lists[index].append( tri_ptr );
00762   tri_sort_list( tri_ptr, index );
00763 
00764   return CUBIT_SUCCESS;
00765 }
00766 
00767 //-------------------------------------------------------------------------
00768 // Function:    insert_at_circumcenter
00769 // Description: insert a node at the circumcenter of a tri
00770 // Author:      chynes
00771 // Date:        6/3/2002
00772 //-------------------------------------------------------------------------
00773 template<class SURF, class TRI, class EDGE, class NODE, class TRICHILD, class NODECHILD, class SIZEFUNC> MY_INLINE
00774 CubitStatus
00775 FacetorUtil<SURF, TRI, EDGE, NODE, TRICHILD, NODECHILD, SIZEFUNC>::insert_at_circumcenter(
00776                                                 TRI* tri_ptr,
00777                                                 DLIList<TRI*>& facet_list,
00778                                                 TRI*& start_tri,
00779                                                 int& curr_visit_flag,
00780                                                 SURF* surface_ptr,
00781                                                 DLIList<TRI*>* sorted_lists,
00782                                                 const int num_lists,
00783                                                 const double interval,
00784                                                 const double quality_angle,
00785                                                 SIZEFUNC* sizing_function,
00786                                                 ParamTool* p_tool)
00787 {
00788   CubitStatus rv = CUBIT_SUCCESS;
00789 
00790   // find the cicumcenter of the triangle and the target size there
00791     //if nodes are colinear do not try to find circumenter
00792   if(!are_nodes_colinear(tri_ptr))
00793   {
00794     PRINT_ERROR("Can't evaluate center of circumcircle\n");
00795     return CUBIT_FAILURE;
00796   }
00797   CubitVector cc = circumcenter( tri_ptr );
00798   
00799   //translate cc into 3D space
00800   CubitVector cc_xyz;
00801   p_tool->transform_to_xyz(cc_xyz, cc);
00802   // get target length in 3D space
00803   double target_length = sizing_function->size_at_point( cc_xyz );
00804   // get new size
00805   double size = get_size(cc, tri_ptr);
00806   // update size
00807   cc.z(size);
00808   // update target_length
00809   target_length = target_length*size;
00810 
00811   // Determine if we should now insert the point.  Allow insertions
00812   // in the general case down to circumcircle size of ALPHA times the
00813   // interpolated target edge length size.  For tris with small
00814   // angles, allow additional inserts to improve the quality down
00815   // to 1/2 the target size
00816   if(!are_nodes_colinear(tri_ptr))
00817   {
00818     PRINT_ERROR("Can't evaluate radius of circumcircle\n");
00819     return CUBIT_FAILURE;
00820   }
00821      
00822     
00823   double r2 = radius( tri_ptr );
00824   CubitBoolean insert = CUBIT_FALSE;
00825   int tsindex = tri_sort_list( tri_ptr );
00826   assert(tsindex > -1);
00827   if (tsindex == 0) 
00828   {
00829     if (r2 > SQR(ALPHA*target_length)) 
00830     {
00831       insert = CUBIT_TRUE;
00832     }
00833   }
00834   else 
00835   {
00836     if (r2 > SQR(0.5*ALPHA*target_length)) 
00837     {
00838       insert = CUBIT_TRUE;
00839     }
00840   }
00841   if (insert) 
00842   {
00843 
00844       // Determine the tris that will be affected by the insertion
00845 
00846     start_tri = tri_ptr;
00847     DLIList <TRI *> tri_list;
00848     rv = natural_neighbor_tris( cc, facet_list,
00849                                                start_tri, surface_ptr,
00850                                                curr_visit_flag, tri_list );
00851       // If it was outside, then we are done with it
00852 
00853     if (tri_list.size() == 0) 
00854     {
00855       return CUBIT_SUCCESS;
00856     }
00857     if (rv != CUBIT_SUCCESS) {
00858       return rv;
00859     }
00860 
00861       // See if we are too close to a boundary
00862 
00863     double x0, y0, x1, y1, cx, cy, edge_radius, dist;
00864     TRI *nntri_ptr;
00865     EDGE *edge_ptr;
00866     int ii, iedge;
00867     for (ii=0; ii<tri_list.size(); ii++) 
00868     {
00869       nntri_ptr = tri_list.get_and_step();
00870       for (iedge=0; iedge<3; iedge++) 
00871       {
00872         edge_ptr = tri_ptr->edge( iedge );
00873 
00874           // An edge encroaches if the distance from the prospective
00875           // new point to the midpoint of the edge is less than 
00876           // half the length of the edge
00877 
00878         if (edge_ptr->marked()) // on the boundary?
00879         {
00880           x0 = (edge_ptr->start_node())->coordinates().x();
00881           y0 = (edge_ptr->start_node())->coordinates().y();
00882           x1 = (edge_ptr->end_node())->coordinates().x();
00883           y1 = (edge_ptr->end_node())->coordinates().y();
00884           cx = (x0 + x1) * 0.5;
00885           cy = (y0 + y1) * 0.5;
00886           edge_radius = sqrt(SQR(x1-x0) + SQR(y1-y0)) * 0.5;       
00887           dist = sqrt( SQR(cx-cc.x()) + SQR(cy-cc.y()) );
00888 
00889             // Edge encroaches: don't insert, return now
00890 
00891           if (dist - edge_radius < FT_INSIDE_TOL) 
00892           {
00893             return CUBIT_SUCCESS;
00894           }
00895         }
00896       }
00897     }
00898 
00899       // Before inserting, remove all the tris on the neighbor
00900       // tri_list from the lists 
00901    
00902     int index;
00903     for (ii=0; ii<tri_list.size(); ii++) 
00904     {
00905       nntri_ptr = tri_list.get_and_step();
00906       index = tri_sort_list( nntri_ptr );
00907       assert(index >= 0);
00908       sorted_lists[index].remove( nntri_ptr );
00909     }
00910 
00911       // Create the new node
00912 
00913     NODE *new_node_ptr = (NODE *)new NODECHILD( cc, surface_ptr );
00914 
00915       // Insert the new node into the mesh
00916 
00917     rv = bowyer_watson_insert( new_node_ptr, tri_list,
00918                                facet_list, curr_visit_flag,
00919                                surface_ptr, start_tri);
00920     if (rv != CUBIT_SUCCESS)
00921         return rv;
00922 
00923       // get the new tris at the node and classify them
00924 
00925     tri_list.clean_out();
00926     new_node_ptr->tris( tri_list );
00927     for (ii=0; ii<tri_list.size() && rv == CUBIT_SUCCESS; ii++) 
00928     {
00929       tri_ptr = tri_list.get_and_step();
00930       rv = classify_tri_by_angle( tri_ptr, sorted_lists, num_lists, interval, quality_angle );
00931     }
00932   }
00933 
00934   return rv;
00935 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines