MOAB: Mesh Oriented datABase  (version 5.3.0)
SmoothFace.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <algorithm>
00004 #include <iomanip>
00005 #include <cassert>
00006 #include <limits>
00007 #include "moab/OrientedBoxTreeTool.hpp"
00008 #include "SmoothFace.hpp"
00009 
00010 #define GEOMETRY_RESABS 1.e-6
00011 #define mbsqr( a )      ( ( a ) * ( a ) )
00012 #define mbcube( a )     ( mbsqr( a ) * ( a ) )
00013 #define mbquart( a )    ( mbsqr( a ) * mbsqr( a ) )
00014 
00015 namespace moab
00016 {
00017 
00018 bool within_tolerance( CartVect& p1, CartVect& p2, const double& tolerance )
00019 {
00020     if( ( fabs( p1[0] - p2[0] ) < tolerance ) && ( fabs( p1[1] - p2[1] ) < tolerance ) &&
00021         ( fabs( p1[2] - p2[2] ) < tolerance ) )
00022         return true;
00023     return false;
00024 }
00025 int numAdjTriInSet( Interface* mb, EntityHandle startEdge, EntityHandle set )
00026 {
00027     std::vector< EntityHandle > adjTri;
00028     mb->get_adjacencies( &startEdge, 1, 2, false, adjTri, Interface::UNION );
00029     // count how many are in the set
00030     int nbInSet = 0;
00031     for( size_t i = 0; i < adjTri.size(); i++ )
00032     {
00033         EntityHandle tri = adjTri[i];
00034         if( mb->contains_entities( set, &tri, 1 ) ) nbInSet++;
00035     }
00036     return nbInSet;
00037 }
00038 
00039 bool debug_surf_eval1 = false;
00040 
00041 SmoothFace::SmoothFace( Interface* mb, EntityHandle surface_set, GeomTopoTool* gTool )
00042     : _markTag( 0 ), _gradientTag( 0 ), _tangentsTag( 0 ), _edgeCtrlTag( 0 ), _facetCtrlTag( 0 ),
00043       _facetEdgeCtrlTag( 0 ), _planeTag( 0 ), _mb( mb ), _set( surface_set ), _my_geomTopoTool( gTool ), _obb_root( 0 ),
00044       _evaluationsCounter( 0 )
00045 {
00046     //_smooth_face = NULL;
00047     //_mbOut->create_meshset(MESHSET_SET, _oSet); //will contain the
00048     // get also the obb_root
00049     if( _my_geomTopoTool )
00050     {
00051         _my_geomTopoTool->get_root( this->_set, _obb_root );
00052         if( debug_surf_eval1 ) _my_geomTopoTool->obb_tree()->stats( _obb_root, std::cout );
00053     }
00054 }
00055 
00056 SmoothFace::~SmoothFace() {}
00057 
00058 double SmoothFace::area()
00059 {
00060     // find the area of this entity
00061     // assert(_smooth_face);
00062     // double area1 = _smooth_face->area();
00063     double totArea = 0.;
00064     for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
00065     {
00066         EntityHandle tria = *it;
00067         const EntityHandle* conn3;
00068         int nnodes;
00069         _mb->get_connectivity( tria, conn3, nnodes );
00070         //
00071         // double coords[9]; // store the coordinates for the nodes
00072         //_mb->get_coords(conn3, 3, coords);
00073         CartVect p[3];
00074         _mb->get_coords( conn3, 3, (double*)&p[0] );
00075         // need to compute the angles
00076         // compute angles and the normal
00077         // CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
00078 
00079         CartVect AB( p[1] - p[0] );  //(p2 - p1);
00080         CartVect BC( p[2] - p[1] );  //(p3 - p2);
00081         CartVect normal = AB * BC;
00082         totArea += normal.length() / 2;
00083     }
00084     return totArea;
00085 }
00086 // these tags will be collected for deletion
00087 void SmoothFace::append_smooth_tags( std::vector< Tag >& smoothTags )
00088 {
00089     // these are created locally, for each smooth face
00090     smoothTags.push_back( _gradientTag );
00091     smoothTags.push_back( _planeTag );
00092 }
00093 void SmoothFace::bounding_box( double box_min[3], double box_max[3] )
00094 {
00095 
00096     for( int i = 0; i < 3; i++ )
00097     {
00098         box_min[i] = _minim[i];
00099         box_max[i] = _maxim[i];
00100     }
00101     // _minim, _maxim
00102 }
00103 
00104 void SmoothFace::move_to_surface( double& x, double& y, double& z )
00105 {
00106 
00107     CartVect loc2( x, y, z );
00108     bool trim    = false;  // is it needed?
00109     bool outside = true;
00110     CartVect closestPoint;
00111 
00112     ErrorCode rval = project_to_facets_main( loc2, trim, outside, &closestPoint, NULL );
00113     if( MB_SUCCESS != rval ) return;
00114     assert( rval == MB_SUCCESS );
00115     x = closestPoint[0];
00116     y = closestPoint[1];
00117     z = closestPoint[2];
00118 }
00119 /*
00120  void SmoothFace::move_to_surface(double& x, double& y, double& z,
00121  double& u_guess, double& v_guess) {
00122  if (debug_surf_eval1) {
00123  std::cout << "move_to_surface called." << std::endl;
00124  }
00125  }*/
00126 
00127 bool SmoothFace::normal_at( double x, double y, double z, double& nx, double& ny, double& nz )
00128 {
00129 
00130     CartVect loc2( x, y, z );
00131 
00132     bool trim    = false;  // is it needed?
00133     bool outside = true;
00134     // CartVect closestPoint;// not needed
00135     // not interested in normal
00136     CartVect normal;
00137     ErrorCode rval = project_to_facets_main( loc2, trim, outside, NULL, &normal );
00138     if( MB_SUCCESS != rval ) return false;
00139     assert( rval == MB_SUCCESS );
00140     nx = normal[0];
00141     ny = normal[1];
00142     nz = normal[2];
00143 
00144     return true;
00145 }
00146 
00147 ErrorCode SmoothFace::compute_control_points_on_edges( double min_dot, Tag edgeCtrlTag, Tag markTag )
00148 {
00149 
00150     _edgeCtrlTag = edgeCtrlTag;
00151     _markTag     = markTag;
00152 
00153     // now, compute control points for all edges that are not marked already (they are no on the
00154     // boundary!)
00155     for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
00156     {
00157         EntityHandle edg = *it;
00158         // is the edge marked? already computed
00159         unsigned char tagVal = 0;
00160         _mb->tag_get_data( _markTag, &edg, 1, &tagVal );
00161         if( tagVal ) continue;
00162         // double min_dot;
00163         init_bezier_edge( edg, min_dot );
00164         tagVal = 1;
00165         _mb->tag_set_data( _markTag, &edg, 1, &tagVal );
00166     }
00167     return MB_SUCCESS;
00168 }
00169 
00170 int SmoothFace::init_gradient()
00171 {
00172     // first, create a Tag for gradient (or normal)
00173     // loop over all triangles in set, and modify the normals according to the angle as weight
00174     // get all the edges from this subset
00175     if( NULL == _mb ) return 1;  // fail
00176     _triangles.clear();
00177     ErrorCode rval = _mb->get_entities_by_type( _set, MBTRI, _triangles );
00178     if( MB_SUCCESS != rval ) return 1;
00179     // get a new range of edges, and decide the loops from here
00180     _edges.clear();
00181     rval = _mb->get_adjacencies( _triangles, 1, true, _edges, Interface::UNION );
00182     assert( rval == MB_SUCCESS );
00183 
00184     rval = _mb->get_adjacencies( _triangles, 0, false, _nodes, Interface::UNION );
00185     assert( rval == MB_SUCCESS );
00186 
00187     // initialize bounding box limits
00188     CartVect vert1;
00189     EntityHandle v1 = _nodes[0];  // first vertex
00190     _mb->get_coords( &v1, 1, (double*)&vert1 );
00191     _minim = vert1;
00192     _maxim = vert1;
00193 
00194     double defNormal[] = { 0., 0., 0. };
00195     // look for a tag name here that is definitely unique. We do not want the tags to interfere with
00196     // each other this normal will be for each node in a face some nodes have multiple normals, if
00197     // they are at the feature edges
00198     unsigned long setId = _mb->id_from_handle( _set );
00199     char name[50]       = { 0 };
00200     sprintf( name, "GRADIENT%lu",
00201              setId );  // name should be something like GRADIENT29, where 29 is the set ID of the face
00202     rval = _mb->tag_get_handle( name, 3, MB_TYPE_DOUBLE, _gradientTag, MB_TAG_DENSE | MB_TAG_CREAT, &defNormal );
00203     assert( rval == MB_SUCCESS );
00204 
00205     double defPlane[4] = { 0., 0., 1., 0. };
00206     // also define a plane tag ; this will be for each triangle
00207     char namePlaneTag[50] = { 0 };
00208     sprintf( namePlaneTag, "PLANE%lu", setId );
00209     rval = _mb->tag_get_handle( "PLANE", 4, MB_TYPE_DOUBLE, _planeTag, MB_TAG_DENSE | MB_TAG_CREAT, &defPlane );
00210     assert( rval == MB_SUCCESS );
00211     // the fourth double is for weight, accumulated at each vertex so far
00212     // maybe not needed in the end
00213     for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
00214     {
00215         EntityHandle tria = *it;
00216         const EntityHandle* conn3;
00217         int nnodes;
00218         _mb->get_connectivity( tria, conn3, nnodes );
00219         if( nnodes != 3 ) return 1;  // error
00220         // double coords[9]; // store the coordinates for the nodes
00221         //_mb->get_coords(conn3, 3, coords);
00222         CartVect p[3];
00223         _mb->get_coords( conn3, 3, (double*)&p[0] );
00224         // need to compute the angles
00225         // compute angles and the normal
00226         // CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
00227 
00228         CartVect AB( p[1] - p[0] );  //(p2 - p1);
00229         CartVect BC( p[2] - p[1] );  //(p3 - p2);
00230         CartVect CA( p[0] - p[2] );  //(p1 - p3);
00231         double a[3];
00232         a[1]            = angle( AB, -BC );  // angle at B (p2), etc.
00233         a[2]            = angle( BC, -CA );
00234         a[0]            = angle( CA, -AB );
00235         CartVect normal = -AB * CA;
00236         normal.normalize();
00237         double plane[4];
00238 
00239         const double* coordNormal = normal.array();
00240 
00241         plane[0] = coordNormal[0];
00242         plane[1] = coordNormal[1];
00243         plane[2] = coordNormal[2];
00244         plane[3] = -normal % p[0];  // dot product
00245         // set the plane
00246         rval = _mb->tag_set_data( _planeTag, &tria, 1, plane );
00247         assert( rval == MB_SUCCESS );
00248 
00249         // add to each vertex the tag value the normal multiplied by the angle
00250         double values[9];
00251 
00252         _mb->tag_get_data( _gradientTag, conn3, 3, values );
00253         for( int i = 0; i < 3; i++ )
00254         {
00255             // values[4*i]+=a[i]; // first is the weight, which we do not really need
00256             values[3 * i + 0] += a[i] * coordNormal[0];
00257             values[3 * i + 1] += a[i] * coordNormal[1];
00258             values[3 * i + 2] += a[i] * coordNormal[2];
00259         }
00260 
00261         // reset those values
00262         _mb->tag_set_data( _gradientTag, conn3, 3, values );
00263     }
00264     // normalize the gradients at each node; maybe not needed here?
00265     // no, we will do it, it is important
00266     int numNodes      = _nodes.size();
00267     double* normalVal = new double[3 * numNodes];
00268     _mb->tag_get_data( _gradientTag, _nodes,
00269                        normalVal );  // get all the normal values at the _nodes
00270     for( int i = 0; i < numNodes; i++ )
00271     {
00272         CartVect p1( &normalVal[3 * i] );
00273         p1.normalize();
00274         p1.get( &normalVal[3 * i] );
00275     }
00276 
00277     // reset the normal values after normalization
00278     _mb->tag_set_data( _gradientTag, _nodes, normalVal );
00279     // print the loops size and some other stuff
00280     if( debug_surf_eval1 )
00281     {
00282         std::cout << " normals at  " << numNodes << " nodes" << std::endl;
00283         int i = 0;
00284         for( Range::iterator it = _nodes.begin(); it != _nodes.end(); ++it, i++ )
00285         {
00286             EntityHandle node = *it;
00287             std::cout << " Node id " << _mb->id_from_handle( node ) << "  " << normalVal[3 * i] << " "
00288                       << normalVal[3 * i + 1] << " " << normalVal[3 * i + 2] << std::endl;
00289         }
00290     }
00291 
00292     delete[] normalVal;
00293 
00294     return 0;
00295 }
00296 
00297 // init bezier edges
00298 // start copy
00299 //===========================================================================
00300 // Function Name: init_bezier_edge
00301 //
00302 // Member Type:  PRIVATE
00303 // Description:  compute the control points for an edge
00304 //===========================================================================
00305 ErrorCode SmoothFace::init_bezier_edge( EntityHandle edge, double )
00306 {
00307     // min dot was used for angle here
00308     // int stat = 0; // CUBIT_SUCCESS;
00309     // all boundaries will be simple, initially
00310     // we may complicate them afterwards
00311 
00312     CartVect ctrl_pts[3];
00313     int nnodes                = 0;
00314     const EntityHandle* conn2 = NULL;
00315     ErrorCode rval            = _mb->get_connectivity( edge, conn2, nnodes );
00316     assert( rval == MB_SUCCESS );
00317     if( MB_SUCCESS != rval ) return rval;
00318 
00319     assert( 2 == nnodes );
00320     // double coords[6]; // store the coordinates for the nodes
00321     CartVect P[2];
00322     // ErrorCode rval = _mb->get_coords(conn2, 2, coords);
00323     rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
00324     assert( rval == MB_SUCCESS );
00325     if( MB_SUCCESS != rval ) return rval;
00326 
00327     // CartVect P0(&coords[0]);
00328     // CartVect P3(&coords[3]);
00329 
00330     // double normalVec[6];
00331     CartVect N[2];
00332     //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
00333     rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
00334     assert( rval == MB_SUCCESS );
00335     if( MB_SUCCESS != rval ) return rval;
00336 
00337     CartVect T[2];  // T0, T3
00338 
00339     rval = _mb->tag_get_data( _tangentsTag, &edge, 1, &T[0] );
00340     assert( rval == MB_SUCCESS );
00341     if( MB_SUCCESS != rval ) return rval;
00342 
00343     rval = init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], ctrl_pts );
00344     assert( rval == MB_SUCCESS );
00345     if( MB_SUCCESS != rval ) return rval;
00346 
00347     rval = _mb->tag_set_data( _edgeCtrlTag, &edge, 1, &ctrl_pts[0] );
00348     assert( rval == MB_SUCCESS );
00349     if( MB_SUCCESS != rval ) return rval;
00350 
00351     if( debug_surf_eval1 )
00352     {
00353         std::cout << "edge: " << _mb->id_from_handle( edge ) << " tangents: " << T[0] << T[1] << std::endl;
00354         std::cout << "  points: " << P[0] << " " << P[1] << std::endl;
00355         std::cout << "  normals: " << N[0] << " " << N[1] << std::endl;
00356         std::cout << "  Control points  " << ctrl_pts[0] << " " << ctrl_pts[1] << " " << ctrl_pts[2] << std::endl;
00357     }
00358 
00359     return MB_SUCCESS;
00360 }
00361 
00362 ErrorCode SmoothFace::compute_tangents_for_each_edge()
00363 // they will be used for control points
00364 {
00365     double defTangents[6] = { 0., 0., 0., 0., 0., 0. };
00366     ErrorCode rval =
00367         _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, _tangentsTag, MB_TAG_DENSE | MB_TAG_CREAT, &defTangents );
00368     if( MB_SUCCESS != rval ) return MB_FAILURE;
00369 
00370     // now, compute Tangents for all edges that are not on boundary, so they are not marked
00371     for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
00372     {
00373         EntityHandle edg = *it;
00374 
00375         int nnodes;
00376         const EntityHandle* conn2;  //
00377         _mb->get_connectivity( edg, conn2, nnodes );
00378         assert( nnodes == 2 );
00379         CartVect P[2];  // store the coordinates for the nodes
00380         rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
00381         if( MB_SUCCESS != rval ) return rval;
00382         assert( rval == MB_SUCCESS );
00383         CartVect T[2];
00384         T[0] = P[1] - P[0];
00385         T[0].normalize();
00386         T[1] = T[0];  //
00387         _mb->tag_set_data( _tangentsTag, &edg, 1,
00388                            (double*)&T[0] );  // set the tangents computed at every edge
00389     }
00390     return MB_SUCCESS;
00391 }
00392 // start copy
00393 //===========================================================================
00394 // Function Name: init_edge_control_points
00395 //
00396 // Member Type:  PRIVATE
00397 // Description:  compute the control points for an edge
00398 //===========================================================================
00399 ErrorCode SmoothFace::init_edge_control_points( CartVect& P0, CartVect& P3, CartVect& N0, CartVect& N3, CartVect& T0,
00400                                                 CartVect& T3, CartVect* Pi )
00401 {
00402     CartVect Vi[4];
00403     Vi[0] = P0;
00404     Vi[3] = P3;
00405     CartVect P03( P3 - P0 );
00406     double di    = P03.length();
00407     double ai    = N0 % N3;  // this is the dot operator, the same as in cgm for CubitVector
00408     double ai0   = N0 % T0;
00409     double ai3   = N3 % T3;
00410     double denom = 4 - ai * ai;
00411     if( fabs( denom ) < 1e-20 )
00412     {
00413         return MB_FAILURE;  // CUBIT_FAILURE;
00414     }
00415     double row   = 6.0e0 * ( 2.0e0 * ai0 + ai * ai3 ) / denom;
00416     double omega = 6.0e0 * ( 2.0e0 * ai3 + ai * ai0 ) / denom;
00417     Vi[1]        = Vi[0] + ( di * ( ( ( 6.0e0 * T0 ) - ( ( 2.0e0 * row ) * N0 ) + ( omega * N3 ) ) / 18.0e0 ) );
00418     Vi[2]        = Vi[3] - ( di * ( ( ( 6.0e0 * T3 ) + ( row * N0 ) - ( ( 2.0e0 * omega ) * N3 ) ) / 18.0e0 ) );
00419     // CartVect Wi[3];
00420     // Wi[0] = Vi[1] - Vi[0];
00421     // Wi[1] = Vi[2] - Vi[1];
00422     // Wi[2] = Vi[3] - Vi[2];
00423 
00424     Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
00425     Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
00426     Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
00427 
00428     return MB_SUCCESS;
00429 }
00430 
00431 ErrorCode SmoothFace::find_edges_orientations( EntityHandle edges[3], const EntityHandle* conn3,
00432                                                int orient[3] )  // maybe we will set it?
00433 {
00434     // find the edge that is adjacent to 2 vertices at a time
00435     for( int i = 0; i < 3; i++ )
00436     {
00437         // edge 0 is 1-2, 1 is 3-1, 2 is 0-1
00438         EntityHandle v[2];
00439         v[0] = conn3[( i + 1 ) % 3];
00440         v[1] = conn3[( i + 2 ) % 3];
00441         std::vector< EntityHandle > adjacencies;
00442         // generate all edges for these two hexes
00443         ErrorCode rval = _mb->get_adjacencies( v, 2, 1, false, adjacencies, Interface::INTERSECT );
00444         if( MB_SUCCESS != rval ) return rval;
00445 
00446         // find the edge connected to both vertices, and then see its orientation
00447         assert( adjacencies.size() == 1 );
00448         const EntityHandle* conn2 = NULL;
00449         int nnodes                = 0;
00450         rval                      = _mb->get_connectivity( adjacencies[0], conn2, nnodes );
00451         assert( rval == MB_SUCCESS );
00452         assert( 2 == nnodes );
00453         edges[i] = adjacencies[0];
00454         // what is the story morning glory?
00455         if( conn2[0] == v[0] && conn2[1] == v[1] )
00456             orient[i] = 1;
00457         else if( conn2[0] == v[1] && conn2[1] == v[0] )
00458             orient[i] = -1;
00459         else
00460             return MB_FAILURE;
00461     }
00462     return MB_SUCCESS;
00463 }
00464 ErrorCode SmoothFace::compute_internal_control_points_on_facets( double, Tag facetCtrlTag, Tag facetEdgeCtrlTag )
00465 {
00466     // collect from each triangle the control points in order
00467     //
00468 
00469     _facetCtrlTag     = facetCtrlTag;
00470     _facetEdgeCtrlTag = facetEdgeCtrlTag;
00471 
00472     for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
00473     {
00474         EntityHandle tri = *it;
00475         // first get connectivity, and the edges
00476         // we need a fast method to retrieve the adjacent edges to each triangle
00477         const EntityHandle* conn3;
00478         int nnodes;
00479         ErrorCode rval = _mb->get_connectivity( tri, conn3, nnodes );
00480         assert( MB_SUCCESS == rval );
00481         if( MB_SUCCESS != rval ) return rval;
00482         assert( 3 == nnodes );
00483 
00484         // would it be easier to do
00485         CartVect vNode[3];  // position at nodes
00486         rval = _mb->get_coords( conn3, 3, (double*)&vNode[0] );
00487         assert( MB_SUCCESS == rval );
00488         if( MB_SUCCESS != rval ) return rval;
00489 
00490         // get gradients (normal) at each node of triangle
00491         CartVect NN[3];
00492         rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
00493         assert( MB_SUCCESS == rval );
00494         if( MB_SUCCESS != rval ) return rval;
00495 
00496         EntityHandle edges[3];
00497         int orient[3];  // + 1 or -1, if the edge is positive or negative within the face
00498         rval = find_edges_orientations( edges, conn3, orient );  // maybe we will set it?
00499         assert( MB_SUCCESS == rval );
00500         if( MB_SUCCESS != rval ) return rval;
00501         // maybe we will store some tags with edges and their orientation with respect to
00502         // a triangle;
00503         CartVect P[3][5];
00504         CartVect N[6], G[6];
00505         // create the linear array for control points on edges, for storage (expensive !!!)
00506         CartVect CP[9];
00507         int index = 0;
00508         // maybe store a tag / entity handle for edges?
00509         for( int i = 0; i < 3; i++ )
00510         {
00511             // populate P and N with the right vectors
00512             int i1       = ( i + 1 ) % 3;  // the first node of the edge
00513             int i2       = ( i + 2 ) % 3;  // the second node of the edge
00514             N[2 * i]     = NN[i1];
00515             N[2 * i + 1] = NN[i2];
00516             P[i][0]      = vNode[i1];
00517             rval         = _mb->tag_get_data( _edgeCtrlTag, &edges[i], 1, &( P[i][1] ) );
00518             // if sense is -1, swap 1 and 3 control points
00519             if( orient[i] == -1 )
00520             {
00521                 CartVect tmp;
00522                 tmp     = P[i][1];
00523                 P[i][1] = P[i][3];
00524                 P[i][3] = tmp;
00525             }
00526             P[i][4] = vNode[i2];
00527             for( int j = 1; j < 4; j++ )
00528                 CP[index++] = P[i][j];
00529 
00530             // the first edge control points
00531         }
00532 
00533         //  stat = facet->get_edge_control_points( P );
00534         init_facet_control_points( N, P, G );
00535         // what do we need to store in the tag control points?
00536         rval = _mb->tag_set_data( _facetCtrlTag, &tri, 1, &G[0] );
00537         assert( MB_SUCCESS == rval );
00538         if( MB_SUCCESS != rval ) return rval;
00539 
00540         // store here again the 9 control points on the edges
00541         rval = _mb->tag_set_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
00542         assert( MB_SUCCESS == rval );
00543         if( MB_SUCCESS != rval ) return rval;
00544         // look at what we retrieve later
00545 
00546         // adjust the bounding box
00547         int j = 0;
00548         for( j = 0; j < 3; j++ )
00549             adjust_bounding_box( vNode[j] );
00550         // edge control points
00551         for( j = 0; j < 9; j++ )
00552             adjust_bounding_box( CP[j] );
00553         // internal facet control points
00554         for( j = 0; j < 6; j++ )
00555             adjust_bounding_box( G[j] );
00556     }
00557     return MB_SUCCESS;
00558 }
00559 void SmoothFace::adjust_bounding_box( CartVect& vect )
00560 {
00561     // _minim, _maxim
00562     for( int j = 0; j < 3; j++ )
00563     {
00564         if( _minim[j] > vect[j] ) _minim[j] = vect[j];
00565         if( _maxim[j] < vect[j] ) _maxim[j] = vect[j];
00566     }
00567 }
00568 //===============================================================
00569 ////Function Name: init_facet_control_points
00570 ////
00571 ////Member Type:  PRIVATE
00572 ////Description:  compute the control points for a facet
00573 ////===============================================================
00574 ErrorCode SmoothFace::init_facet_control_points( CartVect N[6],     // vertex normals (per edge)
00575                                                  CartVect P[3][5],  // edge control points
00576                                                  CartVect G[6] )    // return internal control points
00577 {
00578     CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
00579     double denom;
00580     double lambda[2], mu[2];
00581 
00582     ErrorCode rval = MB_SUCCESS;
00583 
00584     for( int i = 0; i < 3; i++ )
00585     {
00586         N0    = N[i * 2];
00587         N3    = N[i * 2 + 1];
00588         Vi[0] = P[i][0];
00589         Vi[1] = ( P[i][1] - 0.25 * P[i][0] ) / 0.75;
00590         Vi[2] = ( P[i][3] - 0.25 * P[i][4] ) / 0.75;
00591         Vi[3] = P[i][4];
00592         Wi[0] = Vi[1] - Vi[0];
00593         Wi[1] = Vi[2] - Vi[1];
00594         Wi[2] = Vi[3] - Vi[2];
00595         Di[0] = P[( i + 2 ) % 3][3] - 0.5 * ( P[i][1] + P[i][0] );
00596         Di[3] = P[( i + 1 ) % 3][1] - 0.5 * ( P[i][4] + P[i][3] );
00597         Ai[0] = ( N0 * Wi[0] ) / Wi[0].length();
00598         Ai[2] = ( N3 * Wi[2] ) / Wi[2].length();
00599         Ai[1] = Ai[0] + Ai[2];
00600         denom = Ai[1].length();
00601         Ai[1] /= denom;
00602         lambda[0] = ( Di[0] % Wi[0] ) / ( Wi[0] % Wi[0] );
00603         lambda[1] = ( Di[3] % Wi[2] ) / ( Wi[2] % Wi[2] );
00604         mu[0]     = ( Di[0] % Ai[0] );
00605         mu[1]     = ( Di[3] % Ai[2] );
00606         G[i * 2]  = 0.5 * ( P[i][1] + P[i][2] ) + 0.66666666666666 * lambda[0] * Wi[1] +
00607                    0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0] * Ai[1] +
00608                    0.33333333333333 * mu[1] * Ai[0];
00609         G[i * 2 + 1] = 0.5 * ( P[i][2] + P[i][3] ) + 0.33333333333333 * lambda[0] * Wi[2] +
00610                        0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333 * mu[0] * Ai[2] +
00611                        0.66666666666666 * mu[1] * Ai[1];
00612     }
00613     return rval;
00614 }
00615 
00616 void SmoothFace::DumpModelControlPoints()
00617 {
00618     // here, we will dump all control points from edges and facets (6 control points for each facet)
00619     // we may also create some edges; maybe later...
00620     // create a point3D file
00621     // output a Point3D file (special visit format)
00622     unsigned long setId = _mb->id_from_handle( _set );
00623     char name[50]       = { 0 };
00624     sprintf( name, "%lucontrol.Point3D", setId );  // name should be something 2control.Point3D
00625     std::ofstream point3DFile;
00626     point3DFile.open( name );  //("control.Point3D");
00627     point3DFile << "# x y z \n";
00628     std::ofstream point3DEdgeFile;
00629     sprintf( name, "%lucontrolEdge.Point3D", setId );  //
00630     point3DEdgeFile.open( name );                      //("controlEdge.Point3D");
00631     point3DEdgeFile << "# x y z \n";
00632     std::ofstream smoothPoints;
00633     sprintf( name, "%lusmooth.Point3D", setId );  //
00634     smoothPoints.open( name );                    //("smooth.Point3D");
00635     smoothPoints << "# x y z \n";
00636     CartVect controlPoints[3];  // edge control points
00637     for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
00638     {
00639         EntityHandle edge = *it;
00640 
00641         _mb->tag_get_data( _edgeCtrlTag, &edge, 1, (double*)&controlPoints[0] );
00642         for( int i = 0; i < 3; i++ )
00643         {
00644             CartVect& c = controlPoints[i];
00645             point3DEdgeFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
00646         }
00647     }
00648     CartVect controlTriPoints[6];  // triangle control points
00649     CartVect P_facet[3];           // result in 3 "mid" control points
00650     for( Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); ++it2 )
00651     {
00652         EntityHandle tri = *it2;
00653 
00654         _mb->tag_get_data( _facetCtrlTag, &tri, 1, (double*)&controlTriPoints[0] );
00655 
00656         // draw a line of points between pairs of control points
00657         int numPoints = 7;
00658         for( int n = 0; n < numPoints; n++ )
00659         {
00660             double a = 1. * n / ( numPoints - 1 );
00661             double b = 1.0 - a;
00662 
00663             P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
00664             // 1,2,1
00665             P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
00666             // 1,1,2
00667             P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
00668             for( int i = 0; i < 3; i++ )
00669             {
00670                 CartVect& c = P_facet[i];
00671                 point3DFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
00672             }
00673         }
00674 
00675         // evaluate for each triangle a lattice of points
00676         int N = 40;
00677         for( int k = 0; k <= N; k++ )
00678         {
00679             for( int m = 0; m <= N - k; m++ )
00680             {
00681                 int n = N - m - k;
00682                 CartVect areacoord( 1. * k / N, 1. * m / N, 1. * n / N );
00683                 CartVect pt;
00684                 eval_bezier_patch( tri, areacoord, pt );
00685                 smoothPoints << std::setprecision( 11 ) << pt[0] << " " << pt[1] << " " << pt[2] << " \n";
00686             }
00687         }
00688     }
00689     point3DFile.close();
00690     smoothPoints.close();
00691     point3DEdgeFile.close();
00692     return;
00693 }
00694 //===========================================================================
00695 // Function Name: evaluate_single
00696 //
00697 // Member Type:  PUBLIC
00698 // Description:  evaluate edge not associated with a facet (this is used
00699 // by camal edge mesher!!!)
00700 // Note:         t is a value from 0 to 1, for us
00701 //===========================================================================
00702 ErrorCode SmoothFace::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv )
00703 {
00704     CartVect P[2];              // P0 and P1
00705     CartVect controlPoints[3];  // edge control points
00706     double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
00707 
00708     // project the position to the linear edge
00709     // t is from 0 to 1 only!!
00710     // double tt = (t + 1) * 0.5;
00711     if( tt <= 0.0 ) tt = 0.0;
00712     if( tt >= 1.0 ) tt = 1.0;
00713 
00714     int nnodes                = 0;
00715     const EntityHandle* conn2 = NULL;
00716     ErrorCode rval            = _mb->get_connectivity( eh, conn2, nnodes );
00717     assert( rval == MB_SUCCESS );
00718     if( MB_SUCCESS != rval ) return rval;
00719 
00720     rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
00721     assert( rval == MB_SUCCESS );
00722     if( MB_SUCCESS != rval ) return rval;
00723 
00724     rval = _mb->tag_get_data( _edgeCtrlTag, &eh, 1, (double*)&controlPoints[0] );
00725     assert( rval == MB_SUCCESS );
00726     if( MB_SUCCESS != rval ) return rval;
00727 
00728     t2           = tt * tt;
00729     t3           = t2 * tt;
00730     t4           = t3 * tt;
00731     one_minus_t  = 1. - tt;
00732     one_minus_t2 = one_minus_t * one_minus_t;
00733     one_minus_t3 = one_minus_t2 * one_minus_t;
00734     one_minus_t4 = one_minus_t3 * one_minus_t;
00735 
00736     outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
00737            4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
00738 
00739     return MB_SUCCESS;
00740 }
00741 
00742 ErrorCode SmoothFace::eval_bezier_patch( EntityHandle tri, CartVect& areacoord, CartVect& pt )
00743 {
00744     //
00745     // interpolate internal control points
00746 
00747     CartVect gctrl_pts[6];
00748     // get the control points  facet->get_control_points( gctrl_pts );
00749     // init_facet_control_points( N, P, G) ;
00750     // what do we need to store in the tag control points?
00751     ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &tri, 1, &gctrl_pts[0] );  // get all 6 control points
00752     assert( MB_SUCCESS == rval );
00753     if( MB_SUCCESS != rval ) return rval;
00754     const EntityHandle* conn3 = NULL;
00755     int nnodes                = 0;
00756     rval                      = _mb->get_connectivity( tri, conn3, nnodes );
00757     assert( MB_SUCCESS == rval );
00758 
00759     CartVect vN[3];
00760     _mb->get_coords( conn3, 3, (double*)&vN[0] );  // fill the coordinates of the vertices
00761 
00762     if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
00763     {
00764         pt = vN[0];
00765         return MB_SUCCESS;
00766     }
00767     if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
00768     {
00769         pt = vN[0];
00770         return MB_SUCCESS;
00771     }
00772     if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
00773     {
00774         pt = vN[0];
00775         return MB_SUCCESS;
00776     }
00777 
00778     CartVect P_facet[3];
00779     // 2,1,1
00780     P_facet[0] =
00781         ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
00782     // 1,2,1
00783     P_facet[1] =
00784         ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
00785     // 1,1,2
00786     P_facet[2] =
00787         ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
00788 
00789     // sum the contribution from each of the control points
00790 
00791     pt = CartVect( 0. );  // set all to 0, we start adding / accumulating different parts
00792     // first edge is from node 0 to 1, index 2 in
00793 
00794     // retrieve the points, in order, and the control points on edges
00795 
00796     // store here again the 9 control points on the edges
00797     CartVect CP[9];
00798     rval = _mb->tag_get_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
00799     assert( MB_SUCCESS == rval );
00800 
00801     // CubitFacetEdge *edge;
00802     // edge = facet->edge(2);! start with edge 2, from 0-1
00803     int k = 0;
00804     CartVect ctrl_pts[5];
00805     // edge->control_points(facet, ctrl_pts);
00806     ctrl_pts[0] = vN[0];  //
00807     for( k = 1; k < 4; k++ )
00808         ctrl_pts[k] = CP[k + 5];  // for edge index 2
00809     ctrl_pts[4] = vN[1];          //
00810 
00811     // i=4; j=0; k=0;
00812     double B = mbquart( areacoord[0] );
00813     pt += B * ctrl_pts[0];
00814 
00815     // i=3; j=1; k=0;
00816     B = 4.0 * mbcube( areacoord[0] ) * areacoord[1];
00817     pt += B * ctrl_pts[1];
00818 
00819     // i=2; j=2; k=0;
00820     B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[1] );
00821     pt += B * ctrl_pts[2];
00822 
00823     // i=1; j=3; k=0;
00824     B = 4.0 * areacoord[0] * mbcube( areacoord[1] );
00825     pt += B * ctrl_pts[3];
00826 
00827     // edge = facet->edge(0);
00828     // edge->control_points(facet, ctrl_pts);
00829     // edge index 0, from 1 to 2
00830     ctrl_pts[0] = vN[1];  //
00831     for( k = 1; k < 4; k++ )
00832         ctrl_pts[k] = CP[k - 1];  // for edge index 0
00833     ctrl_pts[4] = vN[2];          //
00834 
00835     // i=0; j=4; k=0;
00836     B = mbquart( areacoord[1] );
00837     pt += B * ctrl_pts[0];
00838 
00839     // i=0; j=3; k=1;
00840     B = 4.0 * mbcube( areacoord[1] ) * areacoord[2];
00841     pt += B * ctrl_pts[1];
00842 
00843     // i=0; j=2; k=2;
00844     B = 6.0 * mbsqr( areacoord[1] ) * mbsqr( areacoord[2] );
00845     pt += B * ctrl_pts[2];
00846 
00847     // i=0; j=1; k=3;
00848     B = 4.0 * areacoord[1] * mbcube( areacoord[2] );
00849     pt += B * ctrl_pts[3];
00850 
00851     // edge = facet->edge(1);
00852     // edge->control_points(facet, ctrl_pts);
00853     // edge index 1, from 2 to 0
00854     ctrl_pts[0] = vN[2];  //
00855     for( k = 1; k < 4; k++ )
00856         ctrl_pts[k] = CP[k + 2];  // for edge index 0
00857     ctrl_pts[4] = vN[0];          //
00858 
00859     // i=0; j=0; k=4;
00860     B = mbquart( areacoord[2] );
00861     pt += B * ctrl_pts[0];
00862 
00863     // i=1; j=0; k=3;
00864     B = 4.0 * areacoord[0] * mbcube( areacoord[2] );
00865     pt += B * ctrl_pts[1];
00866 
00867     // i=2; j=0; k=2;
00868     B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[2] );
00869     pt += B * ctrl_pts[2];
00870 
00871     // i=3; j=0; k=1;
00872     B = 4.0 * mbcube( areacoord[0] ) * areacoord[2];
00873     pt += B * ctrl_pts[3];
00874 
00875     // i=2; j=1; k=1;
00876     B = 12.0 * mbsqr( areacoord[0] ) * areacoord[1] * areacoord[2];
00877     pt += B * P_facet[0];
00878 
00879     // i=1; j=2; k=1;
00880     B = 12.0 * areacoord[0] * mbsqr( areacoord[1] ) * areacoord[2];
00881     pt += B * P_facet[1];
00882 
00883     // i=1; j=1; k=2;
00884     B = 12.0 * areacoord[0] * areacoord[1] * mbsqr( areacoord[2] );
00885     pt += B * P_facet[2];
00886 
00887     return MB_SUCCESS;
00888 }
00889 
00890 //===========================================================================
00891 // Function Name: project_to_facet_plane
00892 //
00893 // Member Type:  PUBLIC
00894 // Descriptoin:  Project a point to the plane of a facet
00895 //===========================================================================
00896 void SmoothFace::project_to_facet_plane( EntityHandle tri, CartVect& pt, CartVect& point_on_plane,
00897                                          double& dist_to_plane )
00898 {
00899     double plane[4];
00900 
00901     ErrorCode rval = _mb->tag_get_data( _planeTag, &tri, 1, plane );
00902     if( MB_SUCCESS != rval ) return;
00903     assert( rval == MB_SUCCESS );
00904     // _planeTag
00905     CartVect normal( &plane[0] );  // just first 3 components are used
00906 
00907     double dist    = normal % pt + plane[3];  // coeff d is saved!!!
00908     dist_to_plane  = fabs( dist );
00909     point_on_plane = pt - dist * normal;
00910 
00911     return;
00912 }
00913 
00914 //===========================================================================
00915 // Function Name: facet_area_coordinate
00916 //
00917 // Member Type:  PUBLIC
00918 // Descriptoin:  Determine the area coordinates of a point on the plane
00919 //              of a facet
00920 //===========================================================================
00921 void SmoothFace::facet_area_coordinate( EntityHandle facet, CartVect& pt_on_plane, CartVect& areacoord )
00922 {
00923 
00924     const EntityHandle* conn3 = NULL;
00925     int nnodes                = 0;
00926     ErrorCode rval            = _mb->get_connectivity( facet, conn3, nnodes );
00927     assert( MB_SUCCESS == rval );
00928     if( rval ) {}  // empty statement to prevent compiler warning
00929 
00930     // double coords[9]; // store the coordinates for the nodes
00931     //_mb->get_coords(conn3, 3, coords);
00932     CartVect p[3];
00933     rval = _mb->get_coords( conn3, 3, (double*)&p[0] );
00934     assert( MB_SUCCESS == rval );
00935     if( rval ) {}  // empty statement to prevent compiler warning
00936     double plane[4];
00937 
00938     rval = _mb->tag_get_data( _planeTag, &facet, 1, plane );
00939     assert( rval == MB_SUCCESS );
00940     if( rval ) {}                  // empty statement to prevent compiler warning
00941     CartVect normal( &plane[0] );  // just first 3 components are used
00942 
00943     double area2;
00944 
00945     double tol = GEOMETRY_RESABS * 1.e-5;  // 1.e-11;
00946 
00947     CartVect v1( p[1] - p[0] );
00948     CartVect v2( p[2] - p[0] );
00949 
00950     area2 = ( v1 * v2 ).length_squared();  // the same for CartVect
00951     if( area2 < 100 * tol ) { tol = .01 * area2; }
00952     CartVect absnorm( fabs( normal[0] ), fabs( normal[1] ), fabs( normal[2] ) );
00953 
00954     // project to the closest coordinate plane so we only have to do this in 2D
00955 
00956     if( absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2] )
00957     {
00958         area2 = determ3( p[0][1], p[0][2], p[1][1], p[1][2], p[2][1], p[2][2] );
00959         if( fabs( area2 ) < tol )
00960         {
00961             areacoord = CartVect( -std::numeric_limits< double >::min() );  // .set(
00962                                                                             // -std::numeric_limits<double>::min(),
00963                                                                             // -std::numeric_limits<double>::min(),
00964                                                                             // -std::numeric_limits<double>::min() );
00965         }
00966         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
00967         {
00968             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
00969         }
00970         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
00971         {
00972             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
00973         }
00974         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
00975         {
00976             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
00977         }
00978         else
00979         {
00980 
00981             areacoord[0] = determ3( pt_on_plane[1], pt_on_plane[2], p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
00982 
00983             areacoord[1] = determ3( p[0][1], p[0][2], pt_on_plane[1], pt_on_plane[2], p[2][1], p[2][2] ) / area2;
00984 
00985             areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2], pt_on_plane[1], pt_on_plane[2] ) / area2;
00986         }
00987     }
00988     else if( absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2] )
00989     {
00990 
00991         area2 = determ3( p[0][0], p[0][2], p[1][0], p[1][2], p[2][0], p[2][2] );
00992         if( fabs( area2 ) < tol )
00993         {
00994             areacoord = CartVect( -std::numeric_limits< double >::min() );  //.set(
00995                                                                             //-std::numeric_limits<double>::min(),
00996                                                                             //-std::numeric_limits<double>::min(),
00997                                                                             //-std::numeric_limits<double>::min() );
00998         }
00999         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
01000         {
01001             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
01002         }
01003         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
01004         {
01005             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
01006         }
01007         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
01008         {
01009             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
01010         }
01011         else
01012         {
01013 
01014             areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[2], p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
01015 
01016             areacoord[1] = determ3( p[0][0], p[0][2], pt_on_plane[0], pt_on_plane[2], p[2][0], p[2][2] ) / area2;
01017 
01018             areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2], pt_on_plane[0], pt_on_plane[2] ) / area2;
01019         }
01020     }
01021     else
01022     {
01023         /*area2 = determ3(pt0->x(), pt0->y(),
01024          pt1->x(), pt1->y(),
01025          pt2->x(), pt2->y());*/
01026         area2 = determ3( p[0][0], p[0][1], p[1][0], p[1][1], p[2][0], p[2][1] );
01027         if( fabs( area2 ) < tol )
01028         {
01029             areacoord = CartVect( -std::numeric_limits< double >::min() );  //.set(
01030                                                                             //-std::numeric_limits<double>::min(),
01031                                                                             //-std::numeric_limits<double>::min(),
01032                                                                             //-std::numeric_limits<double>::min() );
01033         }
01034         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
01035         {
01036             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
01037         }
01038         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
01039         {
01040             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
01041         }
01042         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
01043         {
01044             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
01045         }
01046         else
01047         {
01048 
01049             areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[1], p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
01050 
01051             areacoord[1] = determ3( p[0][0], p[0][1], pt_on_plane[0], pt_on_plane[1], p[2][0], p[2][1] ) / area2;
01052 
01053             areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1], pt_on_plane[0], pt_on_plane[1] ) / area2;
01054         }
01055     }
01056 }
01057 
01058 ErrorCode SmoothFace::project_to_facets_main( CartVect& this_point, bool trim, bool& outside,
01059                                               CartVect* closest_point_ptr, CartVect* normal_ptr )
01060 {
01061 
01062     // if there are a lot of facets on this surface - use the OBB search first
01063     // to narrow the selection
01064 
01065     _evaluationsCounter++;
01066     double tolerance = 1.e-5;
01067     std::vector< EntityHandle > facets_out;
01068     // we will start with a list of facets anyway, the best among them wins
01069     ErrorCode rval =
01070         _my_geomTopoTool->obb_tree()->closest_to_location( (double*)&this_point, _obb_root, tolerance, facets_out );
01071     if( MB_SUCCESS != rval ) return rval;
01072 
01073     int interpOrder        = 4;
01074     double compareTol      = 1.e-5;
01075     EntityHandle lastFacet = facets_out.front();
01076     rval = project_to_facets( facets_out, lastFacet, interpOrder, compareTol, this_point, trim, outside,
01077                               closest_point_ptr, normal_ptr );
01078 
01079     return rval;
01080 }
01081 ErrorCode SmoothFace::project_to_facets( std::vector< EntityHandle >& facet_list, EntityHandle& lastFacet,
01082                                          int interpOrder, double compareTol, CartVect& this_point, bool, bool& outside,
01083                                          CartVect* closest_point_ptr, CartVect* normal_ptr )
01084 {
01085 
01086     bool outside_facet      = false;
01087     bool best_outside_facet = true;
01088     double mindist          = 1.e20;
01089     CartVect close_point, best_point( mindist, mindist, mindist ), best_areacoord;
01090     EntityHandle best_facet = 0L;  // no best facet found yet
01091     EntityHandle facet;
01092     assert( facet_list.size() > 0 );
01093 
01094     double big_dist = compareTol * 1.0e3;
01095 
01096     // from the list of close facets, determine the closest point
01097     for( size_t i = 0; i < facet_list.size(); i++ )
01098     {
01099         facet = facet_list[i];
01100         CartVect pt_on_plane;
01101         double dist_to_plane;
01102         project_to_facet_plane( facet, this_point, pt_on_plane, dist_to_plane );
01103 
01104         CartVect areacoord;
01105         // CartVect close_point;
01106         facet_area_coordinate( facet, pt_on_plane, areacoord );
01107         if( interpOrder != 0 )
01108         {
01109 
01110             // modify the areacoord - project to the bezier patch- snaps to the
01111             // edge of the patch if necessary
01112 
01113             if( project_to_facet( facet, this_point, areacoord, close_point, outside_facet, compareTol ) != MB_SUCCESS )
01114             {
01115                 return MB_FAILURE;
01116             }
01117             // if (closest_point_ptr)
01118             //*closest_point_ptr = close_point;
01119         }
01120         // keep track of the minimum distance
01121 
01122         double dist = ( close_point - this_point ).length();  // close_point.distance_between(this_point);
01123         if( ( best_outside_facet == outside_facet && dist < mindist ) ||
01124             ( best_outside_facet && !outside_facet && ( dist < big_dist || best_facet == 0L /*!best_facet*/ ) ) )
01125         {
01126             mindist            = dist;
01127             best_point         = close_point;
01128             best_facet         = facet;
01129             best_areacoord     = areacoord;
01130             best_outside_facet = outside_facet;
01131 
01132             if( dist < compareTol ) { break; }
01133             big_dist = 10.0 * mindist;
01134         }
01135         // facet->marked(1);
01136         // used_facet_list.append(facet);
01137     }
01138 
01139     if( normal_ptr )
01140     {
01141         CartVect normal;
01142         if( eval_bezier_patch_normal( best_facet, best_areacoord, normal ) != MB_SUCCESS ) { return MB_FAILURE; }
01143         *normal_ptr = normal;
01144     }
01145 
01146     if( closest_point_ptr ) { *closest_point_ptr = best_point; }
01147 
01148     outside   = best_outside_facet;
01149     lastFacet = best_facet;
01150 
01151     return MB_SUCCESS;
01152     // end copy
01153 }
01154 
01155 //===========================================================================
01156 // Function Name: project_to_patch
01157 //
01158 // Member Type:  PUBLIC
01159 // Description:  Project a point to a bezier patch. Pass in the areacoord
01160 //              of the point projected to the linear facet.  Function
01161 //              assumes that the point is contained within the patch -
01162 //              if not, it will project to one of its edges.
01163 //===========================================================================
01164 ErrorCode SmoothFace::project_to_patch( EntityHandle facet,   // (IN) the facet where the patch is defined
01165                                         CartVect& ac,         // (IN) area coordinate initial guess (from linear facet)
01166                                         CartVect& pt,         // (IN) point we are projecting to patch
01167                                         CartVect& eval_pt,    // (OUT) The projected point
01168                                         CartVect* eval_norm,  // (OUT) normal at evaluated point
01169                                         bool& outside,        // (OUT) the closest point on patch to pt is on an edge
01170                                         double compare_tol,   // (IN) comparison tolerance
01171                                         int edge_id )         // (IN) only used if this is to be projected to one
01172 //      of the edges.  Otherwise, should be -1
01173 {
01174     ErrorCode status = MB_SUCCESS;
01175 
01176     // see if we are at a vertex
01177 
01178 #define INCR 0.01
01179     const double tol = compare_tol;
01180 
01181     if( is_at_vertex( facet, pt, ac, compare_tol, eval_pt, eval_norm ) )
01182     {
01183         outside = false;
01184         return MB_SUCCESS;
01185     }
01186 
01187     // check if the start ac is inside the patch -if not, then move it there
01188 
01189     int nout          = 0;
01190     const double atol = 0.001;
01191     if( move_ac_inside( ac, atol ) ) nout++;
01192 
01193     int diverge = 0;
01194     int iter    = 0;
01195     CartVect newpt;
01196     eval_bezier_patch( facet, ac, newpt );
01197     CartVect move   = pt - newpt;
01198     double lastdist = move.length();
01199     double bestdist = lastdist;
01200     CartVect bestac = ac;
01201     CartVect bestpt = newpt;
01202     CartVect bestnorm( 0, 0, 0 );
01203 
01204     // If we are already close enough, then return now
01205 
01206     if( lastdist <= tol && !eval_norm && nout == 0 )
01207     {
01208         eval_pt = pt;
01209         outside = false;
01210         return status;
01211     }
01212 
01213     double ratio, mag, umove, vmove, det, distnew, movedist;
01214     CartVect lastpt = newpt;
01215     CartVect lastac = ac;
01216     CartVect norm;
01217     CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
01218     CartVect du, dv, newac;
01219     bool done = false;
01220     while( !done )
01221     {
01222 
01223         // We will be locating the projected point within the u,v,w coordinate
01224         // system of the triangular bezier patch.  Since u+v+w=1, the components
01225         // are not linearly independent.  We will choose only two of the
01226         // coordinates to use and compute the third.
01227 
01228         int system;
01229         if( lastac[0] >= lastac[1] && lastac[0] >= lastac[2] ) { system = 0; }
01230         else if( lastac[1] >= lastac[2] )
01231         {
01232             system = 1;
01233         }
01234         else
01235         {
01236             system = 2;
01237         }
01238 
01239         // compute the surface derivatives with respect to each
01240         // of the barycentric coordinates
01241 
01242         if( system == 1 || system == 2 )
01243         {
01244             xac[0] = lastac[0] + INCR;  // xac.x( lastac.x() + INCR );
01245             if( lastac[1] + lastac[2] == 0.0 ) return MB_FAILURE;
01246             ratio  = lastac[2] / ( lastac[1] + lastac[2] );  // ratio = lastac.z() / (lastac.y() + lastac.z());
01247             xac[1] = ( 1.0 - xac[0] ) * ( 1.0 - ratio );     // xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
01248             xac[2] = 1.0 - xac[0] - xac[1];                  // xac.z( 1.0 - xac.x() - xac.y() );
01249             eval_bezier_patch( facet, xac, xpt );
01250             xvec = xpt - lastpt;
01251             xvec /= INCR;
01252         }
01253         if( system == 0 || system == 2 )
01254         {
01255             yac[1] = ( lastac[1] + INCR );      // yac.y( lastac.y() + INCR );
01256             if( lastac[0] + lastac[2] == 0.0 )  // if (lastac.x() + lastac.z() == 0.0)
01257                 return MB_FAILURE;
01258             ratio  = lastac[2] / ( lastac[0] + lastac[2] );   // ratio = lastac.z() / (lastac.x() + lastac.z());
01259             yac[0] = ( ( 1.0 - yac[1] ) * ( 1.0 - ratio ) );  // yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
01260             yac[2] = ( 1.0 - yac[0] - yac[1] );               // yac.z( 1.0 - yac.x() - yac.y() );
01261             eval_bezier_patch( facet, yac, ypt );
01262             yvec = ypt - lastpt;
01263             yvec /= INCR;
01264         }
01265         if( system == 0 || system == 1 )
01266         {
01267             zac[2] = ( lastac[2] + INCR );      // zac.z( lastac.z() + INCR );
01268             if( lastac[0] + lastac[1] == 0.0 )  // if (lastac.x() + lastac.y() == 0.0)
01269                 return MB_FAILURE;
01270             ratio  = lastac[1] / ( lastac[0] + lastac[1] );   // ratio = lastac.y() / (lastac.x() + lastac.y());
01271             zac[0] = ( ( 1.0 - zac[2] ) * ( 1.0 - ratio ) );  // zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
01272             zac[1] = ( 1.0 - zac[0] - zac[2] );               // zac.y( 1.0 - zac.x() - zac.z() );
01273             eval_bezier_patch( facet, zac, zpt );
01274             zvec = zpt - lastpt;
01275             zvec /= INCR;
01276         }
01277 
01278         // compute the surface normal
01279 
01280         switch( system )
01281         {
01282             case 0:
01283                 du = yvec;
01284                 dv = zvec;
01285                 break;
01286             case 1:
01287                 du = zvec;
01288                 dv = xvec;
01289                 break;
01290             case 2:
01291                 du = xvec;
01292                 dv = yvec;
01293                 break;
01294         }
01295         norm = du * dv;
01296         mag  = norm.length();
01297         if( mag < DBL_EPSILON )
01298         {
01299             return MB_FAILURE;
01300             // do something else here (it is likely a flat triangle -
01301             // so try evaluating just an edge of the bezier patch)
01302         }
01303         norm /= mag;
01304         if( iter == 0 ) bestnorm = norm;
01305 
01306         // project the move vector to the tangent plane
01307 
01308         move = ( norm * move ) * norm;
01309 
01310         // compute an equivalent u-v-w vector
01311 
01312         CartVect absnorm( fabs( norm[0] ), fabs( norm[1] ), fabs( norm[2] ) );
01313         if( absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0] )
01314         {
01315             det = du[0] * dv[1] - dv[0] * du[1];
01316             if( fabs( det ) <= DBL_EPSILON )
01317             {
01318                 return MB_FAILURE;  // do something else here
01319             }
01320             umove = ( move[0] * dv[1] - dv[0] * move[1] ) / det;
01321             vmove = ( du[0] * move[1] - move[0] * du[1] ) / det;
01322         }
01323         else if( absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0] )
01324         {
01325             det = du[0] * dv[2] - dv[0] * du[2];
01326             if( fabs( det ) <= DBL_EPSILON ) { return MB_FAILURE; }
01327             umove = ( move[0] * dv[2] - dv[0] * move[2] ) / det;
01328             vmove = ( du[0] * move[2] - move[0] * du[2] ) / det;
01329         }
01330         else
01331         {
01332             det = du[1] * dv[2] - dv[1] * du[2];
01333             if( fabs( det ) <= DBL_EPSILON ) { return MB_FAILURE; }
01334             umove = ( move[1] * dv[2] - dv[1] * move[2] ) / det;
01335             vmove = ( du[1] * move[2] - move[1] * du[2] ) / det;
01336         }
01337 
01338         /* === compute the new u-v coords and evaluate surface at new location */
01339 
01340         switch( system )
01341         {
01342             case 0:
01343                 newac[1] = ( lastac[1] + umove );          // newac.y( lastac.y() + umove );
01344                 newac[2] = ( lastac[2] + vmove );          // newac.z( lastac.z() + vmove );
01345                 newac[0] = ( 1.0 - newac[1] - newac[2] );  // newac.x( 1.0 - newac.y() - newac.z() );
01346                 break;
01347             case 1:
01348                 newac[2] = ( lastac[2] + umove );          // newac.z( lastac.z() + umove );
01349                 newac[0] = ( lastac[0] + vmove );          // newac.x( lastac.x() + vmove );
01350                 newac[1] = ( 1.0 - newac[2] - newac[0] );  // newac.y( 1.0 - newac.z() - newac.x() );
01351                 break;
01352             case 2:
01353                 newac[0] = ( lastac[0] + umove );          // newac.x( lastac.x() + umove );
01354                 newac[1] = ( lastac[1] + vmove );          // newac.y( lastac.y() + vmove );
01355                 newac[2] = ( 1.0 - newac[0] - newac[1] );  // newac.z( 1.0 - newac.x() - newac.y() );
01356                 break;
01357         }
01358 
01359         // Keep it inside the patch
01360 
01361         if( newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol ) { nout = 0; }
01362         else
01363         {
01364             if( move_ac_inside( newac, atol ) ) nout++;
01365         }
01366 
01367         // Evaluate at the new location
01368 
01369         if( edge_id != -1 ) ac_at_edge( newac, newac, edge_id );  // move to edge first
01370         eval_bezier_patch( facet, newac, newpt );
01371 
01372         // Check for convergence
01373 
01374         distnew  = ( pt - newpt ).length();  // pt.distance_between(newpt);
01375         move     = newpt - lastpt;
01376         movedist = move.length();
01377         if( movedist < tol || distnew < tol )
01378         {
01379             done = true;
01380             if( distnew < bestdist )
01381             {
01382                 bestdist = distnew;
01383                 bestac   = newac;
01384                 bestpt   = newpt;
01385                 bestnorm = norm;
01386             }
01387         }
01388         else
01389         {
01390 
01391             // don't allow more than 30 iterations
01392 
01393             iter++;
01394             if( iter > 30 )
01395             {
01396                 // if (movedist > tol * 100.0) nout=1;
01397                 done = true;
01398             }
01399 
01400             // Check for divergence - don't allow more than 5 divergent
01401             // iterations
01402 
01403             if( distnew > lastdist )
01404             {
01405                 diverge++;
01406                 if( diverge > 10 )
01407                 {
01408                     done = true;
01409                     // if (movedist > tol * 100.0) nout=1;
01410                 }
01411             }
01412 
01413             // Check if we are continuing to project outside the facet.
01414             // If so, then stop now
01415 
01416             if( nout > 3 ) { done = true; }
01417 
01418             // set up for next iteration
01419 
01420             if( !done )
01421             {
01422                 if( distnew < bestdist )
01423                 {
01424                     bestdist = distnew;
01425                     bestac   = newac;
01426                     bestpt   = newpt;
01427                     bestnorm = norm;
01428                 }
01429                 lastdist = distnew;
01430                 lastpt   = newpt;
01431                 lastac   = newac;
01432                 move     = pt - lastpt;
01433             }
01434         }
01435     }
01436 
01437     eval_pt = bestpt;
01438     if( eval_norm ) { *eval_norm = bestnorm; }
01439     outside = ( nout > 0 ) ? true : false;
01440     ac      = bestac;
01441 
01442     return status;
01443 }
01444 
01445 //===========================================================================
01446 // Function Name: ac_at_edge
01447 //
01448 // Member Type:  PRIVATE
01449 // Description:  determine the area coordinate of the facet at the edge
01450 //===========================================================================
01451 void SmoothFace::ac_at_edge( CartVect& fac,  // facet area coordinate
01452                              CartVect& eac,  // edge area coordinate
01453                              int edge_id )   // id of edge
01454 {
01455     double u, v, w;
01456     switch( edge_id )
01457     {
01458         case 0:
01459             u = 0.0;
01460             v = fac[1] / ( fac[1] + fac[2] );  // v = fac.y() / (fac.y() + fac.z());
01461             w = 1.0 - v;
01462             break;
01463         case 1:
01464             u = fac[0] / ( fac[0] + fac[2] );  // u = fac.x() / (fac.x() + fac.z());
01465             v = 0.0;
01466             w = 1.0 - u;
01467             break;
01468         case 2:
01469             u = fac[0] / ( fac[0] + fac[1] );  // u = fac.x() / (fac.x() + fac.y());
01470             v = 1.0 - u;
01471             w = 0.0;
01472             break;
01473         default:
01474             assert( 0 );
01475             u = -1;  // needed to eliminate warnings about used before set
01476             v = -1;  // needed to eliminate warnings about used before set
01477             w = -1;  // needed to eliminate warnings about used before set
01478             break;
01479     }
01480     eac[0] = u;
01481     eac[1] = v;
01482     eac[2] = w;  //= CartVect(u, v, w);
01483 }
01484 
01485 //===========================================================================
01486 // Function Name: project_to_facet
01487 //
01488 // Member Type:  PUBLIC
01489 // Description:  project to a single facet.  Uses the input areacoord as
01490 //              a starting guess.
01491 //===========================================================================
01492 ErrorCode SmoothFace::project_to_facet( EntityHandle facet, CartVect& pt, CartVect& areacoord, CartVect& close_point,
01493                                         bool& outside_facet, double compare_tol )
01494 {
01495     const EntityHandle* conn3 = NULL;
01496     int nnodes                = 0;
01497     _mb->get_connectivity( facet, conn3, nnodes );
01498     //
01499     // double coords[9]; // store the coordinates for the nodes
01500     //_mb->get_coords(conn3, 3, coords);
01501     CartVect p[3];
01502     _mb->get_coords( conn3, 3, (double*)&p[0] );
01503 
01504     int edge_id    = -1;
01505     ErrorCode stat = project_to_patch( facet, areacoord, pt, close_point, NULL, outside_facet, compare_tol, edge_id );
01506     /* }
01507      break;
01508      }*/
01509 
01510     return stat;
01511 }
01512 
01513 //===========================================================================
01514 // Function Name: is_at_vertex
01515 //
01516 // Member Type:  PRIVATE
01517 // Description:  determine if the point is at one of the facet's vertices
01518 //===========================================================================
01519 bool SmoothFace::is_at_vertex( EntityHandle facet,        // (IN) facet we are evaluating
01520                                CartVect& pt,              // (IN) the point
01521                                CartVect& ac,              // (IN) the ac of the point on the facet plane
01522                                double compare_tol,        // (IN) return TRUE of closer than this
01523                                CartVect& eval_pt,         // (OUT) location at vertex if TRUE
01524                                CartVect* eval_norm_ptr )  // (OUT) normal at vertex if TRUE
01525 {
01526     double dist;
01527     CartVect vert_loc;
01528     const double actol = 0.1;
01529     // get coordinates get_coords
01530     const EntityHandle* conn3 = NULL;
01531     int nnodes                = 0;
01532     _mb->get_connectivity( facet, conn3, nnodes );
01533     //
01534     // double coords[9]; // store the coordinates for the nodes
01535     //_mb->get_coords(conn3, 3, coords);
01536     CartVect p[3];
01537     _mb->get_coords( conn3, 3, (double*)&p[0] );
01538     // also get the normals at nodes
01539     CartVect NN[3];
01540     _mb->tag_get_data( _gradientTag, conn3, 3, (double*)&NN[0] );
01541     if( fabs( ac[0] ) < actol && fabs( ac[1] ) < actol )
01542     {
01543         vert_loc = p[2];
01544         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01545         if( dist <= compare_tol )
01546         {
01547             eval_pt = vert_loc;
01548             if( eval_norm_ptr ) { *eval_norm_ptr = NN[2]; }
01549             return true;
01550         }
01551     }
01552 
01553     if( fabs( ac[0] ) < actol && fabs( ac[2] ) < actol )
01554     {
01555         vert_loc = p[1];
01556         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01557         if( dist <= compare_tol )
01558         {
01559             eval_pt = vert_loc;
01560             if( eval_norm_ptr )
01561             {
01562                 *eval_norm_ptr = NN[1];  // facet->point(1)->normal( facet );
01563             }
01564             return true;
01565         }
01566     }
01567 
01568     if( fabs( ac[1] ) < actol && fabs( ac[2] ) < actol )
01569     {
01570         vert_loc = p[0];
01571         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01572         if( dist <= compare_tol )
01573         {
01574             eval_pt = vert_loc;
01575             if( eval_norm_ptr ) { *eval_norm_ptr = NN[0]; }
01576             return true;
01577         }
01578     }
01579 
01580     return false;
01581 }
01582 
01583 //===========================================================================
01584 // Function Name: move_ac_inside
01585 //
01586 // Member Type:  PRIVATE
01587 // Description:  find the closest area coordinate to the boundary of the
01588 //              patch if any of its components are < 0
01589 //              Return if the ac was modified.
01590 //===========================================================================
01591 bool SmoothFace::move_ac_inside( CartVect& ac, double tol )
01592 {
01593     int nout = 0;
01594     if( ac[0] < -tol )
01595     {
01596         ac[0] = 0.0;
01597         ac[1] = ac[1] / ( ac[1] + ac[2] );  //( ac.y() / (ac.y() + ac.z()) ;
01598         ac[2] = 1. - ac[1];                 // ac.z( 1.0 - ac.y() );
01599         nout++;
01600     }
01601     if( ac[1] < -tol )
01602     {
01603         ac[1] = 0.;                         // ac.y( 0.0 );
01604         ac[0] = ac[0] / ( ac[0] + ac[2] );  // ac.x( ac.x() / (ac.x() + ac.z()) );
01605         ac[2] = 1. - ac[0];                 // ac.z( 1.0 - ac.x() );
01606         nout++;
01607     }
01608     if( ac[2] < -tol )
01609     {
01610         ac[2] = 0.;                         // ac.z( 0.0 );
01611         ac[0] = ac[0] / ( ac[0] + ac[1] );  // ac.x( ac.x() / (ac.x() + ac.y()) );
01612         ac[1] = 1. - ac[0];                 // ac.y( 1.0 - ac.x() );
01613         nout++;
01614     }
01615     return ( nout > 0 ) ? true : false;
01616 }
01617 
01618 //===========================================================================
01619 // Function Name: hodograph
01620 //
01621 // Member Type:  PUBLIC
01622 // Description:  get the hodograph control points for the facet
01623 // Note:  This is a triangle cubic patch that is defined by the
01624 //       normals of quartic facet control point lattice.  Returned coordinates
01625 //       in Nijk are defined by the following diagram
01626 //
01627 //
01628 //                         *9               index  polar
01629 //                        / \                 0     300    point(0)
01630 //                       /   \                1     210
01631 //                     7*-----*8              2     120
01632 //                     / \   / \              3     030    point(1)
01633 //                    /   \ /   \             4     201
01634 //                  4*----5*-----*6           5     111
01635 //                  / \   / \   / \           6     021
01636 //                 /   \ /   \ /   \          7     102
01637 //                *-----*-----*-----*         8     012
01638 //                0     1     2     3         9     003    point(2)
01639 //
01640 
01641 //===========================================================================
01642 // Function Name: eval_bezier_patch_normal
01643 //
01644 // Member Type:  PRIVATE
01645 // Description:  evaluate the Bezier patch defined at a facet
01646 //===========================================================================
01647 ErrorCode SmoothFace::eval_bezier_patch_normal( EntityHandle facet, CartVect& areacoord, CartVect& normal )
01648 {
01649     // interpolate internal control points
01650 
01651     CartVect gctrl_pts[6];
01652     // facet->get_control_points( gctrl_pts );
01653     ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &facet, 1, &gctrl_pts[0] );
01654     assert( rval == MB_SUCCESS );
01655     if( MB_SUCCESS != rval ) return rval;
01656     // _gradientTag
01657     // get normals at points
01658     const EntityHandle* conn3 = NULL;
01659     int nnodes                = 0;
01660     rval                      = _mb->get_connectivity( facet, conn3, nnodes );
01661     if( MB_SUCCESS != rval ) return rval;
01662 
01663     CartVect NN[3];
01664     rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
01665     assert( rval == MB_SUCCESS );
01666     if( MB_SUCCESS != rval ) return rval;
01667 
01668     if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
01669     {
01670         normal = NN[0];
01671         return MB_SUCCESS;
01672     }
01673     if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
01674     {
01675         normal = NN[1];  // facet->point(1)->normal(facet);
01676         return MB_SUCCESS;
01677     }
01678     if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
01679     {
01680         normal = NN[2];  // facet->point(2)->normal(facet);
01681         return MB_SUCCESS;
01682     }
01683 
01684     // compute the hodograph of the quartic Gregory patch
01685 
01686     CartVect Nijk[10];
01687     // hodograph(facet,areacoord,Nijk);
01688     // start copy from hodograph
01689     // CubitVector gctrl_pts[6];
01690     // facet->get_control_points( gctrl_pts );
01691     CartVect P_facet[3];
01692 
01693     // 2,1,1
01694     /*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
01695      (areacoord.y() * gctrl_pts[3] +
01696      areacoord.z() * gctrl_pts[4]);*/
01697     P_facet[0] =
01698         ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
01699     // 1,2,1
01700     /*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
01701      (areacoord.x() * gctrl_pts[0] +
01702      areacoord.z() * gctrl_pts[5]);*/
01703     P_facet[1] =
01704         ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
01705     // 1,1,2
01706     /*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
01707      (areacoord.x() * gctrl_pts[1] +
01708      areacoord.y() * gctrl_pts[2]);*/
01709     P_facet[2] =
01710         ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
01711 
01712     // corner control points are just the normals at the points
01713 
01714     // 3, 0, 0
01715     Nijk[0] = NN[0];
01716     // 0, 3, 0
01717     Nijk[3] = NN[1];
01718     // 0, 0, 3
01719     Nijk[9] = NN[2];  // facet->point(2)->normal(facet);
01720 
01721     // fill in the boundary control points.  Define as the normal to the local
01722     // triangle formed by the quartic control point lattice
01723 
01724     // store here again the 9 control points on the edges
01725     CartVect CP[9];  // 9 control points on the edges,
01726     rval = _mb->tag_get_data( _facetEdgeCtrlTag, &facet, 1, &CP[0] );
01727     if( MB_SUCCESS != rval ) return rval;
01728     // there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
01729     // CubitFacetEdge *edge;
01730     // edge = facet->edge( 2 );
01731     // CubitVector ctrl_pts[5];
01732     // edge->control_points(facet, ctrl_pts);
01733 
01734     // 2, 1, 0
01735     // Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
01736     Nijk[1] = ( CP[7] - CP[6] ) * ( P_facet[0] - CP[6] );
01737     Nijk[1].normalize();
01738 
01739     // 1, 2, 0
01740     // Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
01741     Nijk[2] = ( CP[8] - CP[7] ) * ( P_facet[1] - CP[7] );
01742     Nijk[2].normalize();
01743 
01744     // edge = facet->edge( 0 );
01745     // edge->control_points(facet, ctrl_pts);
01746 
01747     // 0, 2, 1
01748     // Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
01749     Nijk[6] = ( CP[0] - P_facet[1] ) * ( CP[1] - P_facet[1] );
01750     Nijk[6].normalize();
01751 
01752     // 0, 1, 2
01753     // Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
01754     Nijk[8] = ( CP[1] - P_facet[2] ) * ( CP[2] - P_facet[2] );
01755     Nijk[8].normalize();
01756 
01757     // edge = facet->edge( 1 );
01758     // edge->control_points(facet, ctrl_pts);
01759 
01760     // 1, 0, 2
01761     // Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
01762     Nijk[7] = ( P_facet[2] - CP[4] ) * ( CP[3] - CP[4] );
01763     Nijk[7].normalize();
01764 
01765     // 2, 0, 1
01766     // Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
01767     Nijk[4] = ( P_facet[0] - CP[5] ) * ( CP[4] - CP[5] );
01768     Nijk[4].normalize();
01769 
01770     // 1, 1, 1
01771     Nijk[5] = ( P_facet[1] - P_facet[0] ) * ( P_facet[2] - P_facet[0] );
01772     Nijk[5].normalize();
01773     // end copy from hodograph
01774 
01775     // sum the contribution from each of the control points
01776 
01777     normal = CartVect( 0.0e0, 0.0e0, 0.0e0 );
01778 
01779     // i=3; j=0; k=0;
01780     // double Bsum = 0.0;
01781     double B = mbcube( areacoord[0] );
01782     // Bsum += B;
01783     normal += B * Nijk[0];
01784 
01785     // i=2; j=1; k=0;
01786     B = 3.0 * mbsqr( areacoord[0] ) * areacoord[1];
01787     // Bsum += B;
01788     normal += B * Nijk[1];
01789 
01790     // i=1; j=2; k=0;
01791     B = 3.0 * areacoord[0] * mbsqr( areacoord[1] );
01792     // Bsum += B;
01793     normal += B * Nijk[2];
01794 
01795     // i=0; j=3; k=0;
01796     B = mbcube( areacoord[1] );
01797     // Bsum += B;
01798     normal += B * Nijk[3];
01799 
01800     // i=2; j=0; k=1;
01801     B = 3.0 * mbsqr( areacoord[0] ) * areacoord[2];
01802     // Bsum += B;
01803     normal += B * Nijk[4];
01804 
01805     // i=1; j=1; k=1;
01806     B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
01807     // Bsum += B;
01808     normal += B * Nijk[5];
01809 
01810     // i=0; j=2; k=1;
01811     B = 3.0 * mbsqr( areacoord[1] ) * areacoord[2];
01812     // Bsum += B;
01813     normal += B * Nijk[6];
01814 
01815     // i=1; j=0; k=2;
01816     B = 3.0 * areacoord[0] * mbsqr( areacoord[2] );
01817     // Bsum += B;
01818     normal += B * Nijk[7];
01819 
01820     // i=0; j=1; k=2;
01821     B = 3.0 * areacoord[1] * mbsqr( areacoord[2] );
01822     // Bsum += B;
01823     normal += B * Nijk[8];
01824 
01825     // i=0; j=0; k=3;
01826     B = mbcube( areacoord[2] );
01827     // Bsum += B;
01828     normal += B * Nijk[9];
01829 
01830     // assert(fabs(Bsum - 1.0) < 1e-9);
01831 
01832     normal.normalize();
01833 
01834     return MB_SUCCESS;
01835 }
01836 
01837 ErrorCode SmoothFace::get_normals_for_vertices( const EntityHandle* conn2, CartVect N[2] )
01838 // this method will be called to retrieve the normals needed in the calculation of control edge
01839 // points..
01840 {
01841     // CartVect N[2];
01842     //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
01843     ErrorCode rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
01844     return rval;
01845 }
01846 
01847 ErrorCode SmoothFace::ray_intersection_correct( EntityHandle,       // (IN) the facet where the patch is defined
01848                                                 CartVect& pt,       // (IN) shoot from
01849                                                 CartVect& ray,      // (IN) ray direction
01850                                                 CartVect& eval_pt,  // (INOUT) The intersection point
01851                                                 double& distance,   // (IN OUT) the new distance
01852                                                 bool& outside )
01853 {
01854     // find a point on the smooth surface
01855     CartVect currentPoint = eval_pt;
01856     int numIter           = 0;
01857     double improvement    = 1.e20;
01858     CartVect diff;
01859     while( numIter++ < 5 && improvement > 0.01 )
01860     {
01861         CartVect newPos;
01862 
01863         bool trim = false;  // is it needed?
01864         outside   = true;
01865         CartVect closestPoint;
01866         CartVect normal;
01867 
01868         ErrorCode rval = project_to_facets_main( currentPoint, trim, outside, &newPos, &normal );
01869         if( MB_SUCCESS != rval ) return rval;
01870         assert( rval == MB_SUCCESS );
01871         diff        = newPos - currentPoint;
01872         improvement = diff.length();
01873         // ( pt + t * ray - closest ) % normal = 0;
01874         // intersect tangent plane that goes through closest point with the direction
01875         // t = normal%(closest-pt) / normal%ray;
01876         double dot = normal % ray;  // if it is 0, get out while we can
01877         if( dot < 0.00001 )
01878         {
01879             // bad convergence, get out, do not modify anything
01880             return MB_SUCCESS;
01881         }
01882         double t     = ( ( newPos - pt ) % normal ) / ( dot );
01883         currentPoint = pt + t * ray;
01884     }
01885     eval_pt  = currentPoint;
01886     diff     = currentPoint - pt;
01887     distance = diff.length();
01888     return MB_SUCCESS;
01889 }
01890 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines