cgma
GeoTet.cpp
Go to the documentation of this file.
00001 //-------------------------------------------------------------------------
00002 // Filename      : GeoTet.cpp
00003 //
00004 // Purpose       : Tet class used for geometry operations.  See also GeoNode.
00005 //
00006 // Description   : light-weight entity used for computational geometry 
00007 //                 tools.  Don't load this up with extra stuff!
00008 //
00009 // Creator       : Steve Owen
00010 //
00011 // Creation Date : 8/13/2003
00012 //
00013 // Owner         : Steve Owen
00014 //-------------------------------------------------------------------------
00015 
00016 #include "GeoTet.hpp"
00017 #include "GeoNode.hpp"
00018 #include "GfxDebug.hpp"
00019 #include "CubitMessage.hpp"
00020 
00021 #ifndef SQR
00022 #define SQR(x) ((x) * (x))
00023 #endif
00024 
00025 //--------------------------------------------------------------------------
00026 // Function: GeoTet
00027 // Description: constructor
00028 // Author: sjowen
00029 // Date: 8/13/2003
00030 //---------------------------------------------------------------------------
00031 GeoTet::GeoTet( GeoNode *nodes[4] ) : isMarked(0), isInside(0)
00032 {
00033   mNodes[0] = nodes[0];
00034   mNodes[1] = nodes[1];
00035   mNodes[2] = nodes[2];
00036   mNodes[3] = nodes[3];
00037 }
00038 
00039 //--------------------------------------------------------------------------
00040 // Function: GeoTet
00041 // Description: destructor
00042 // Author: sjowen
00043 // Date: 8/13/2003
00044 //---------------------------------------------------------------------------
00045 GeoTet::~GeoTet()
00046 {
00047 }
00048 
00049 //--------------------------------------------------------------------------
00050 // Function: tet_nodes
00051 // Description: return the nodes on this tet
00052 // Author: sjowen
00053 // Date: 8/13/2003
00054 //---------------------------------------------------------------------------
00055 void GeoTet::tet_nodes( GeoNode *&na, GeoNode *&nb, GeoNode *&nc, GeoNode *&nd )
00056 {
00057   na = mNodes[0];
00058   nb = mNodes[1];
00059   nc = mNodes[2];
00060   nd = mNodes[3];
00061 }
00062 
00063 //--------------------------------------------------------------------------
00064 // Function: tet_face_nodes
00065 // Description: return the on a face of this tet.
00066 // Author: mlstate
00067 // Date: 9/08/2004
00068 //---------------------------------------------------------------------------
00069 void GeoTet::tet_face_nodes( int face_index, GeoNode *&na, GeoNode *&nb, GeoNode *&nc )
00070 {
00071     int node_face_idx[4][3] = { { 1, 2, 3 },
00072                                 { 3, 2, 0 },
00073                                 { 0, 1, 3 },
00074                                 { 2, 1, 0 } };
00075     na = mNodes[node_face_idx[face_index][0]];
00076     nb = mNodes[node_face_idx[face_index][1]];
00077     nc = mNodes[node_face_idx[face_index][2]];
00078 
00079 }
00080 
00081 //--------------------------------------------------------------------------
00082 // Function: get_connected_tet
00083 // Description: return the adjacent tet at the face of face_indx.
00084 // Author: mlstate
00085 // Date: 9/08/2004
00086 //---------------------------------------------------------------------------
00087 GeoTet *GeoTet::get_connected_tet( int face_indx )
00088 {
00089     GeoNode *n1, *n2, *n3;
00090     tet_face_nodes( face_indx, n1, n2, n3 );
00091     return get_connected_tet( n1, n2, n3 );
00092 }
00093 
00094 //--------------------------------------------------------------------------
00095 // Function: get_connected_tet
00096 // Description: return the adjacent tet at the face defined by the three nodes
00097 // Author: sjowen
00098 // Date: 8/13/2003
00099 //---------------------------------------------------------------------------
00100 GeoTet *GeoTet::get_connected_tet( GeoNode *na, GeoNode *nb, GeoNode *nc )
00101 {
00102   GeoTet *adj_tet = NULL;
00103   int ii;
00104 
00105   // get the tets adjacent 
00106   DLIList<GeoTet *> *tet_list_ptr = na->tet_list();
00107   GeoTet *tet_ptr;
00108   for (ii=0; ii<tet_list_ptr->size() && !adj_tet; ii++)
00109   {
00110     tet_ptr = tet_list_ptr->get_and_step();
00111     if (tet_ptr != this)
00112     {
00113       if (tet_ptr->node_index( nb ) >= 0 &&
00114           tet_ptr->node_index( nc ) >= 0)
00115       {
00116         adj_tet = tet_ptr;
00117       }
00118     }
00119   }
00120 
00121   return adj_tet;
00122 }
00123 
00124 //--------------------------------------------------------------------------
00125 // Function: node_index
00126 // Description: return the index of the node in the tet
00127 // Author: sjowen
00128 // Date: 8/13/2003
00129 //---------------------------------------------------------------------------
00130 int GeoTet::node_index( GeoNode *node_ptr )
00131 {
00132   int ii;
00133   for (ii=0; ii<4; ii++)
00134   {
00135     if (mNodes[ii] == node_ptr)
00136       return ii;
00137   }
00138   return -1;
00139 }
00140 
00141 //--------------------------------------------------------------------------
00142 // Function: opposite_edge
00143 // Description: return nodes on the opposite edge of the tet.  a_node and 
00144 //              b_node must be on this tet.  c_node and d_node are in
00145 //              no particular order
00146 // Author: sjowen
00147 // Date: 01/29/2004
00148 //---------------------------------------------------------------------------
00149 void GeoTet::opposite_edge( GeoNode *a_node, GeoNode *b_node, 
00150                             GeoNode *&c_node, GeoNode *&d_node )
00151 {
00152   int ii;
00153   c_node = d_node = NULL;
00154   for (ii=0; ii<4; ii++)
00155   {
00156     if (mNodes[ii] != a_node && mNodes[ii] != b_node)
00157     {
00158       if (!c_node)
00159       {
00160         c_node = mNodes[ii];
00161       }
00162       else if(!d_node)
00163       {
00164         d_node = mNodes[ii];
00165       }
00166       else
00167       {
00168         PRINT_ERROR("a_node or b_node are not on this tet.\n");
00169         return;
00170       }
00171     }
00172   }
00173 }
00174 
00175 //--------------------------------------------------------------------------
00176 // Function: dgtet
00177 // Description: global debug function to draw a GeoTet 
00178 // Author: sjowen
00179 // Date: 8/13/2003
00180 //---------------------------------------------------------------------------
00181 void dgtet( GeoTet *tet )
00182 {
00183   GfxDebug::draw_geotet( tet, CUBIT_YELLOW_INDEX );
00184   GfxDebug::flush();
00185 }
00186 
00187 //--------------------------------------------------------------------------
00188 // Function: dgnode
00189 // Description: global debug function to draw a GeoNode 
00190 // Author: sjowen
00191 // Date: 8/13/2003
00192 //---------------------------------------------------------------------------
00193 void dgnode( GeoNode *node )
00194 {
00195   GfxDebug::draw_geonode( node, CUBIT_RED_INDEX );
00196   GfxDebug::flush();
00197 }
00198 
00199 CubitStatus GeoTet::circumsphere( CubitVector &center, double *radius )
00200 {
00201   double reltol = DBL_EPSILON * 100.0;
00202 
00203   CubitVector a = mNodes[0]->coordinates();
00204   CubitVector b = mNodes[1]->coordinates();
00205   CubitVector c = mNodes[2]->coordinates();
00206   CubitVector d = mNodes[3]->coordinates();
00207 
00208   CubitVector da = a - d;
00209   CubitVector db = b - d;
00210   CubitVector dc = c - d;
00211 
00212   double rhsa = 0.5*(SQR(da.x()) + SQR(da.y()) + SQR(da.z()));
00213   double rhsb = 0.5*(SQR(db.x()) + SQR(db.y()) + SQR(db.z()));
00214   double rhsc = 0.5*(SQR(dc.x()) + SQR(dc.y()) + SQR(dc.z()));
00215 
00216   double cpa = db.y()*dc.z() - dc.y()*db.z();
00217   double cpb = dc.y()*da.z() - da.y()*dc.z();
00218   double cpc = da.y()*db.z() - db.y()*da.z();
00219 
00220   double det = da.x()*cpa + db.x()*cpb + dc.x()*cpc;
00221 
00222   double xmax = CUBIT_MAX(fabs(a.x()),fabs(b.x()));
00223          xmax = CUBIT_MAX(xmax,fabs(c.x())); 
00224          xmax = CUBIT_MAX(xmax,fabs(d.x()));
00225   double ymax = CUBIT_MAX(fabs(a.y()),fabs(b.y()));
00226          ymax = CUBIT_MAX(ymax,fabs(c.y())); 
00227          ymax = CUBIT_MAX(ymax,fabs(d.y()));
00228   double zmax = CUBIT_MAX(fabs(a.z()),fabs(b.z()));
00229          zmax = CUBIT_MAX(zmax,fabs(c.z())); 
00230          zmax = CUBIT_MAX(zmax,fabs(d.z()));
00231   double tolabs = reltol*xmax*ymax*zmax;
00232   if (fabs(det) <= tolabs) {
00233     return CUBIT_FAILURE;
00234   }
00235   center.x( (rhsa*cpa + rhsb*cpb + rhsc*cpc)/det );
00236   cpa = db.x()*rhsc - dc.x()*rhsb;
00237   cpb = dc.x()*rhsa - da.x()*rhsc;
00238   cpc = da.x()*rhsb - db.x()*rhsa;
00239   center.y( (da.z()*cpa + db.z()*cpb + dc.z()*cpc)/det );
00240   center.z( -(da.y()*cpa + db.y()*cpb + dc.y()*cpc)/det );
00241   center += d;
00242 
00243   if ( radius )
00244   {
00245       double radsq = SQR(center.x()) + SQR(center.y()) + SQR(center.z());
00246       *radius = sqrt( radsq );
00247   }
00248 
00249   return CUBIT_SUCCESS;
00250 }
00251 
00252 // EOF
00253 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines