Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00400                                                 CartVect& P3,
00401                                                 CartVect& N0,
00402                                                 CartVect& N3,
00403                                                 CartVect& T0,
00404                                                 CartVect& T3,
00405                                                 CartVect* Pi )
00406 {
00407     CartVect Vi[4];
00408     Vi[0] = P0;
00409     Vi[3] = P3;
00410     CartVect P03( P3 - P0 );
00411     double di    = P03.length();
00412     double ai    = N0 % N3;  // this is the dot operator, the same as in cgm for CubitVector
00413     double ai0   = N0 % T0;
00414     double ai3   = N3 % T3;
00415     double denom = 4 - ai * ai;
00416     if( fabs( denom ) < 1e-20 )
00417     {
00418         return MB_FAILURE;  // CUBIT_FAILURE;
00419     }
00420     double row   = 6.0e0 * ( 2.0e0 * ai0 + ai * ai3 ) / denom;
00421     double omega = 6.0e0 * ( 2.0e0 * ai3 + ai * ai0 ) / denom;
00422     Vi[1]        = Vi[0] + ( di * ( ( ( 6.0e0 * T0 ) - ( ( 2.0e0 * row ) * N0 ) + ( omega * N3 ) ) / 18.0e0 ) );
00423     Vi[2]        = Vi[3] - ( di * ( ( ( 6.0e0 * T3 ) + ( row * N0 ) - ( ( 2.0e0 * omega ) * N3 ) ) / 18.0e0 ) );
00424     // CartVect Wi[3];
00425     // Wi[0] = Vi[1] - Vi[0];
00426     // Wi[1] = Vi[2] - Vi[1];
00427     // Wi[2] = Vi[3] - Vi[2];
00428 
00429     Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
00430     Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
00431     Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
00432 
00433     return MB_SUCCESS;
00434 }
00435 
00436 ErrorCode SmoothFace::find_edges_orientations( EntityHandle edges[3],
00437                                                const EntityHandle* conn3,
00438                                                int orient[3] )  // maybe we will set it?
00439 {
00440     // find the edge that is adjacent to 2 vertices at a time
00441     for( int i = 0; i < 3; i++ )
00442     {
00443         // edge 0 is 1-2, 1 is 3-1, 2 is 0-1
00444         EntityHandle v[2];
00445         v[0] = conn3[( i + 1 ) % 3];
00446         v[1] = conn3[( i + 2 ) % 3];
00447         std::vector< EntityHandle > adjacencies;
00448         // generate all edges for these two hexes
00449         ErrorCode rval = _mb->get_adjacencies( v, 2, 1, false, adjacencies, Interface::INTERSECT );
00450         if( MB_SUCCESS != rval ) return rval;
00451 
00452         // find the edge connected to both vertices, and then see its orientation
00453         assert( adjacencies.size() == 1 );
00454         const EntityHandle* conn2 = NULL;
00455         int nnodes                = 0;
00456         rval                      = _mb->get_connectivity( adjacencies[0], conn2, nnodes );
00457         assert( rval == MB_SUCCESS );
00458         assert( 2 == nnodes );
00459         edges[i] = adjacencies[0];
00460         // what is the story morning glory?
00461         if( conn2[0] == v[0] && conn2[1] == v[1] )
00462             orient[i] = 1;
00463         else if( conn2[0] == v[1] && conn2[1] == v[0] )
00464             orient[i] = -1;
00465         else
00466             return MB_FAILURE;
00467     }
00468     return MB_SUCCESS;
00469 }
00470 ErrorCode SmoothFace::compute_internal_control_points_on_facets( double, Tag facetCtrlTag, Tag facetEdgeCtrlTag )
00471 {
00472     // collect from each triangle the control points in order
00473     //
00474 
00475     _facetCtrlTag     = facetCtrlTag;
00476     _facetEdgeCtrlTag = facetEdgeCtrlTag;
00477 
00478     for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
00479     {
00480         EntityHandle tri = *it;
00481         // first get connectivity, and the edges
00482         // we need a fast method to retrieve the adjacent edges to each triangle
00483         const EntityHandle* conn3;
00484         int nnodes;
00485         ErrorCode rval = _mb->get_connectivity( tri, conn3, nnodes );
00486         assert( MB_SUCCESS == rval );
00487         if( MB_SUCCESS != rval ) return rval;
00488         assert( 3 == nnodes );
00489 
00490         // would it be easier to do
00491         CartVect vNode[3];  // position at nodes
00492         rval = _mb->get_coords( conn3, 3, (double*)&vNode[0] );
00493         assert( MB_SUCCESS == rval );
00494         if( MB_SUCCESS != rval ) return rval;
00495 
00496         // get gradients (normal) at each node of triangle
00497         CartVect NN[3];
00498         rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
00499         assert( MB_SUCCESS == rval );
00500         if( MB_SUCCESS != rval ) return rval;
00501 
00502         EntityHandle edges[3];
00503         int orient[3];  // + 1 or -1, if the edge is positive or negative within the face
00504         rval = find_edges_orientations( edges, conn3, orient );  // maybe we will set it?
00505         assert( MB_SUCCESS == rval );
00506         if( MB_SUCCESS != rval ) return rval;
00507         // maybe we will store some tags with edges and their orientation with respect to
00508         // a triangle;
00509         CartVect P[3][5];
00510         CartVect N[6], G[6];
00511         // create the linear array for control points on edges, for storage (expensive !!!)
00512         CartVect CP[9];
00513         int index = 0;
00514         // maybe store a tag / entity handle for edges?
00515         for( int i = 0; i < 3; i++ )
00516         {
00517             // populate P and N with the right vectors
00518             int i1       = ( i + 1 ) % 3;  // the first node of the edge
00519             int i2       = ( i + 2 ) % 3;  // the second node of the edge
00520             N[2 * i]     = NN[i1];
00521             N[2 * i + 1] = NN[i2];
00522             P[i][0]      = vNode[i1];
00523             rval         = _mb->tag_get_data( _edgeCtrlTag, &edges[i], 1, &( P[i][1] ) );
00524             // if sense is -1, swap 1 and 3 control points
00525             if( orient[i] == -1 )
00526             {
00527                 CartVect tmp;
00528                 tmp     = P[i][1];
00529                 P[i][1] = P[i][3];
00530                 P[i][3] = tmp;
00531             }
00532             P[i][4] = vNode[i2];
00533             for( int j = 1; j < 4; j++ )
00534                 CP[index++] = P[i][j];
00535 
00536             // the first edge control points
00537         }
00538 
00539         //  stat = facet->get_edge_control_points( P );
00540         init_facet_control_points( N, P, G );
00541         // what do we need to store in the tag control points?
00542         rval = _mb->tag_set_data( _facetCtrlTag, &tri, 1, &G[0] );
00543         assert( MB_SUCCESS == rval );
00544         if( MB_SUCCESS != rval ) return rval;
00545 
00546         // store here again the 9 control points on the edges
00547         rval = _mb->tag_set_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
00548         assert( MB_SUCCESS == rval );
00549         if( MB_SUCCESS != rval ) return rval;
00550         // look at what we retrieve later
00551 
00552         // adjust the bounding box
00553         int j = 0;
00554         for( j = 0; j < 3; j++ )
00555             adjust_bounding_box( vNode[j] );
00556         // edge control points
00557         for( j = 0; j < 9; j++ )
00558             adjust_bounding_box( CP[j] );
00559         // internal facet control points
00560         for( j = 0; j < 6; j++ )
00561             adjust_bounding_box( G[j] );
00562     }
00563     return MB_SUCCESS;
00564 }
00565 void SmoothFace::adjust_bounding_box( CartVect& vect )
00566 {
00567     // _minim, _maxim
00568     for( int j = 0; j < 3; j++ )
00569     {
00570         if( _minim[j] > vect[j] ) _minim[j] = vect[j];
00571         if( _maxim[j] < vect[j] ) _maxim[j] = vect[j];
00572     }
00573 }
00574 //===============================================================
00575 ////Function Name: init_facet_control_points
00576 ////
00577 ////Member Type:  PRIVATE
00578 ////Description:  compute the control points for a facet
00579 ////===============================================================
00580 ErrorCode SmoothFace::init_facet_control_points( CartVect N[6],     // vertex normals (per edge)
00581                                                  CartVect P[3][5],  // edge control points
00582                                                  CartVect G[6] )    // return internal control points
00583 {
00584     CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
00585     double denom;
00586     double lambda[2], mu[2];
00587 
00588     ErrorCode rval = MB_SUCCESS;
00589 
00590     for( int i = 0; i < 3; i++ )
00591     {
00592         N0    = N[i * 2];
00593         N3    = N[i * 2 + 1];
00594         Vi[0] = P[i][0];
00595         Vi[1] = ( P[i][1] - 0.25 * P[i][0] ) / 0.75;
00596         Vi[2] = ( P[i][3] - 0.25 * P[i][4] ) / 0.75;
00597         Vi[3] = P[i][4];
00598         Wi[0] = Vi[1] - Vi[0];
00599         Wi[1] = Vi[2] - Vi[1];
00600         Wi[2] = Vi[3] - Vi[2];
00601         Di[0] = P[( i + 2 ) % 3][3] - 0.5 * ( P[i][1] + P[i][0] );
00602         Di[3] = P[( i + 1 ) % 3][1] - 0.5 * ( P[i][4] + P[i][3] );
00603         Ai[0] = ( N0 * Wi[0] ) / Wi[0].length();
00604         Ai[2] = ( N3 * Wi[2] ) / Wi[2].length();
00605         Ai[1] = Ai[0] + Ai[2];
00606         denom = Ai[1].length();
00607         Ai[1] /= denom;
00608         lambda[0] = ( Di[0] % Wi[0] ) / ( Wi[0] % Wi[0] );
00609         lambda[1] = ( Di[3] % Wi[2] ) / ( Wi[2] % Wi[2] );
00610         mu[0]     = ( Di[0] % Ai[0] );
00611         mu[1]     = ( Di[3] % Ai[2] );
00612         G[i * 2]  = 0.5 * ( P[i][1] + P[i][2] ) + 0.66666666666666 * lambda[0] * Wi[1] +
00613                    0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0] * Ai[1] +
00614                    0.33333333333333 * mu[1] * Ai[0];
00615         G[i * 2 + 1] = 0.5 * ( P[i][2] + P[i][3] ) + 0.33333333333333 * lambda[0] * Wi[2] +
00616                        0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333 * mu[0] * Ai[2] +
00617                        0.66666666666666 * mu[1] * Ai[1];
00618     }
00619     return rval;
00620 }
00621 
00622 void SmoothFace::DumpModelControlPoints()
00623 {
00624     // here, we will dump all control points from edges and facets (6 control points for each facet)
00625     // we may also create some edges; maybe later...
00626     // create a point3D file
00627     // output a Point3D file (special visit format)
00628     unsigned long setId = _mb->id_from_handle( _set );
00629     char name[50]       = { 0 };
00630     sprintf( name, "%lucontrol.Point3D", setId );  // name should be something 2control.Point3D
00631     std::ofstream point3DFile;
00632     point3DFile.open( name );  //("control.Point3D");
00633     point3DFile << "# x y z \n";
00634     std::ofstream point3DEdgeFile;
00635     sprintf( name, "%lucontrolEdge.Point3D", setId );  //
00636     point3DEdgeFile.open( name );                      //("controlEdge.Point3D");
00637     point3DEdgeFile << "# x y z \n";
00638     std::ofstream smoothPoints;
00639     sprintf( name, "%lusmooth.Point3D", setId );  //
00640     smoothPoints.open( name );                    //("smooth.Point3D");
00641     smoothPoints << "# x y z \n";
00642     CartVect controlPoints[3];  // edge control points
00643     for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
00644     {
00645         EntityHandle edge = *it;
00646 
00647         _mb->tag_get_data( _edgeCtrlTag, &edge, 1, (double*)&controlPoints[0] );
00648         for( int i = 0; i < 3; i++ )
00649         {
00650             CartVect& c = controlPoints[i];
00651             point3DEdgeFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
00652         }
00653     }
00654     CartVect controlTriPoints[6];  // triangle control points
00655     CartVect P_facet[3];           // result in 3 "mid" control points
00656     for( Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); ++it2 )
00657     {
00658         EntityHandle tri = *it2;
00659 
00660         _mb->tag_get_data( _facetCtrlTag, &tri, 1, (double*)&controlTriPoints[0] );
00661 
00662         // draw a line of points between pairs of control points
00663         int numPoints = 7;
00664         for( int n = 0; n < numPoints; n++ )
00665         {
00666             double a = 1. * n / ( numPoints - 1 );
00667             double b = 1.0 - a;
00668 
00669             P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
00670             // 1,2,1
00671             P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
00672             // 1,1,2
00673             P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
00674             for( int i = 0; i < 3; i++ )
00675             {
00676                 CartVect& c = P_facet[i];
00677                 point3DFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
00678             }
00679         }
00680 
00681         // evaluate for each triangle a lattice of points
00682         int N = 40;
00683         for( int k = 0; k <= N; k++ )
00684         {
00685             for( int m = 0; m <= N - k; m++ )
00686             {
00687                 int n = N - m - k;
00688                 CartVect areacoord( 1. * k / N, 1. * m / N, 1. * n / N );
00689                 CartVect pt;
00690                 eval_bezier_patch( tri, areacoord, pt );
00691                 smoothPoints << std::setprecision( 11 ) << pt[0] << " " << pt[1] << " " << pt[2] << " \n";
00692             }
00693         }
00694     }
00695     point3DFile.close();
00696     smoothPoints.close();
00697     point3DEdgeFile.close();
00698     return;
00699 }
00700 //===========================================================================
00701 // Function Name: evaluate_single
00702 //
00703 // Member Type:  PUBLIC
00704 // Description:  evaluate edge not associated with a facet (this is used
00705 // by camal edge mesher!!!)
00706 // Note:         t is a value from 0 to 1, for us
00707 //===========================================================================
00708 ErrorCode SmoothFace::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv )
00709 {
00710     CartVect P[2];              // P0 and P1
00711     CartVect controlPoints[3];  // edge control points
00712     double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
00713 
00714     // project the position to the linear edge
00715     // t is from 0 to 1 only!!
00716     // double tt = (t + 1) * 0.5;
00717     if( tt <= 0.0 ) tt = 0.0;
00718     if( tt >= 1.0 ) tt = 1.0;
00719 
00720     int nnodes                = 0;
00721     const EntityHandle* conn2 = NULL;
00722     ErrorCode rval            = _mb->get_connectivity( eh, conn2, nnodes );
00723     assert( rval == MB_SUCCESS );
00724     if( MB_SUCCESS != rval ) return rval;
00725 
00726     rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
00727     assert( rval == MB_SUCCESS );
00728     if( MB_SUCCESS != rval ) return rval;
00729 
00730     rval = _mb->tag_get_data( _edgeCtrlTag, &eh, 1, (double*)&controlPoints[0] );
00731     assert( rval == MB_SUCCESS );
00732     if( MB_SUCCESS != rval ) return rval;
00733 
00734     t2           = tt * tt;
00735     t3           = t2 * tt;
00736     t4           = t3 * tt;
00737     one_minus_t  = 1. - tt;
00738     one_minus_t2 = one_minus_t * one_minus_t;
00739     one_minus_t3 = one_minus_t2 * one_minus_t;
00740     one_minus_t4 = one_minus_t3 * one_minus_t;
00741 
00742     outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
00743            4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
00744 
00745     return MB_SUCCESS;
00746 }
00747 
00748 ErrorCode SmoothFace::eval_bezier_patch( EntityHandle tri, CartVect& areacoord, CartVect& pt )
00749 {
00750     //
00751     // interpolate internal control points
00752 
00753     CartVect gctrl_pts[6];
00754     // get the control points  facet->get_control_points( gctrl_pts );
00755     // init_facet_control_points( N, P, G) ;
00756     // what do we need to store in the tag control points?
00757     ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &tri, 1, &gctrl_pts[0] );  // get all 6 control points
00758     assert( MB_SUCCESS == rval );
00759     if( MB_SUCCESS != rval ) return rval;
00760     const EntityHandle* conn3 = NULL;
00761     int nnodes                = 0;
00762     rval                      = _mb->get_connectivity( tri, conn3, nnodes );
00763     assert( MB_SUCCESS == rval );
00764 
00765     CartVect vN[3];
00766     _mb->get_coords( conn3, 3, (double*)&vN[0] );  // fill the coordinates of the vertices
00767 
00768     if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
00769     {
00770         pt = vN[0];
00771         return MB_SUCCESS;
00772     }
00773     if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
00774     {
00775         pt = vN[0];
00776         return MB_SUCCESS;
00777     }
00778     if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
00779     {
00780         pt = vN[0];
00781         return MB_SUCCESS;
00782     }
00783 
00784     CartVect P_facet[3];
00785     // 2,1,1
00786     P_facet[0] =
00787         ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
00788     // 1,2,1
00789     P_facet[1] =
00790         ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
00791     // 1,1,2
00792     P_facet[2] =
00793         ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
00794 
00795     // sum the contribution from each of the control points
00796 
00797     pt = CartVect( 0. );  // set all to 0, we start adding / accumulating different parts
00798     // first edge is from node 0 to 1, index 2 in
00799 
00800     // retrieve the points, in order, and the control points on edges
00801 
00802     // store here again the 9 control points on the edges
00803     CartVect CP[9];
00804     rval = _mb->tag_get_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
00805     assert( MB_SUCCESS == rval );
00806 
00807     // CubitFacetEdge *edge;
00808     // edge = facet->edge(2);! start with edge 2, from 0-1
00809     int k = 0;
00810     CartVect ctrl_pts[5];
00811     // edge->control_points(facet, ctrl_pts);
00812     ctrl_pts[0] = vN[0];  //
00813     for( k = 1; k < 4; k++ )
00814         ctrl_pts[k] = CP[k + 5];  // for edge index 2
00815     ctrl_pts[4] = vN[1];          //
00816 
00817     // i=4; j=0; k=0;
00818     double B = mbquart( areacoord[0] );
00819     pt += B * ctrl_pts[0];
00820 
00821     // i=3; j=1; k=0;
00822     B = 4.0 * mbcube( areacoord[0] ) * areacoord[1];
00823     pt += B * ctrl_pts[1];
00824 
00825     // i=2; j=2; k=0;
00826     B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[1] );
00827     pt += B * ctrl_pts[2];
00828 
00829     // i=1; j=3; k=0;
00830     B = 4.0 * areacoord[0] * mbcube( areacoord[1] );
00831     pt += B * ctrl_pts[3];
00832 
00833     // edge = facet->edge(0);
00834     // edge->control_points(facet, ctrl_pts);
00835     // edge index 0, from 1 to 2
00836     ctrl_pts[0] = vN[1];  //
00837     for( k = 1; k < 4; k++ )
00838         ctrl_pts[k] = CP[k - 1];  // for edge index 0
00839     ctrl_pts[4] = vN[2];          //
00840 
00841     // i=0; j=4; k=0;
00842     B = mbquart( areacoord[1] );
00843     pt += B * ctrl_pts[0];
00844 
00845     // i=0; j=3; k=1;
00846     B = 4.0 * mbcube( areacoord[1] ) * areacoord[2];
00847     pt += B * ctrl_pts[1];
00848 
00849     // i=0; j=2; k=2;
00850     B = 6.0 * mbsqr( areacoord[1] ) * mbsqr( areacoord[2] );
00851     pt += B * ctrl_pts[2];
00852 
00853     // i=0; j=1; k=3;
00854     B = 4.0 * areacoord[1] * mbcube( areacoord[2] );
00855     pt += B * ctrl_pts[3];
00856 
00857     // edge = facet->edge(1);
00858     // edge->control_points(facet, ctrl_pts);
00859     // edge index 1, from 2 to 0
00860     ctrl_pts[0] = vN[2];  //
00861     for( k = 1; k < 4; k++ )
00862         ctrl_pts[k] = CP[k + 2];  // for edge index 0
00863     ctrl_pts[4] = vN[0];          //
00864 
00865     // i=0; j=0; k=4;
00866     B = mbquart( areacoord[2] );
00867     pt += B * ctrl_pts[0];
00868 
00869     // i=1; j=0; k=3;
00870     B = 4.0 * areacoord[0] * mbcube( areacoord[2] );
00871     pt += B * ctrl_pts[1];
00872 
00873     // i=2; j=0; k=2;
00874     B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[2] );
00875     pt += B * ctrl_pts[2];
00876 
00877     // i=3; j=0; k=1;
00878     B = 4.0 * mbcube( areacoord[0] ) * areacoord[2];
00879     pt += B * ctrl_pts[3];
00880 
00881     // i=2; j=1; k=1;
00882     B = 12.0 * mbsqr( areacoord[0] ) * areacoord[1] * areacoord[2];
00883     pt += B * P_facet[0];
00884 
00885     // i=1; j=2; k=1;
00886     B = 12.0 * areacoord[0] * mbsqr( areacoord[1] ) * areacoord[2];
00887     pt += B * P_facet[1];
00888 
00889     // i=1; j=1; k=2;
00890     B = 12.0 * areacoord[0] * areacoord[1] * mbsqr( areacoord[2] );
00891     pt += B * P_facet[2];
00892 
00893     return MB_SUCCESS;
00894 }
00895 
00896 //===========================================================================
00897 // Function Name: project_to_facet_plane
00898 //
00899 // Member Type:  PUBLIC
00900 // Descriptoin:  Project a point to the plane of a facet
00901 //===========================================================================
00902 void SmoothFace::project_to_facet_plane( EntityHandle tri,
00903                                          CartVect& pt,
00904                                          CartVect& point_on_plane,
00905                                          double& dist_to_plane )
00906 {
00907     double plane[4];
00908 
00909     ErrorCode rval = _mb->tag_get_data( _planeTag, &tri, 1, plane );
00910     if( MB_SUCCESS != rval ) return;
00911     assert( rval == MB_SUCCESS );
00912     // _planeTag
00913     CartVect normal( &plane[0] );  // just first 3 components are used
00914 
00915     double dist    = normal % pt + plane[3];  // coeff d is saved!!!
00916     dist_to_plane  = fabs( dist );
00917     point_on_plane = pt - dist * normal;
00918 
00919     return;
00920 }
00921 
00922 //===========================================================================
00923 // Function Name: facet_area_coordinate
00924 //
00925 // Member Type:  PUBLIC
00926 // Descriptoin:  Determine the area coordinates of a point on the plane
00927 //              of a facet
00928 //===========================================================================
00929 void SmoothFace::facet_area_coordinate( EntityHandle facet, CartVect& pt_on_plane, CartVect& areacoord )
00930 {
00931 
00932     const EntityHandle* conn3 = NULL;
00933     int nnodes                = 0;
00934     ErrorCode rval            = _mb->get_connectivity( facet, conn3, nnodes );
00935     assert( MB_SUCCESS == rval );
00936     if( rval )
00937     {
00938     }  // empty statement to prevent compiler warning
00939 
00940     // double coords[9]; // store the coordinates for the nodes
00941     //_mb->get_coords(conn3, 3, coords);
00942     CartVect p[3];
00943     rval = _mb->get_coords( conn3, 3, (double*)&p[0] );
00944     assert( MB_SUCCESS == rval );
00945     if( rval )
00946     {
00947     }  // empty statement to prevent compiler warning
00948     double plane[4];
00949 
00950     rval = _mb->tag_get_data( _planeTag, &facet, 1, plane );
00951     assert( rval == MB_SUCCESS );
00952     if( rval )
00953     {
00954     }                              // empty statement to prevent compiler warning
00955     CartVect normal( &plane[0] );  // just first 3 components are used
00956 
00957     double area2;
00958 
00959     double tol = GEOMETRY_RESABS * 1.e-5;  // 1.e-11;
00960 
00961     CartVect v1( p[1] - p[0] );
00962     CartVect v2( p[2] - p[0] );
00963 
00964     area2 = ( v1 * v2 ).length_squared();  // the same for CartVect
00965     if( area2 < 100 * tol )
00966     {
00967         tol = .01 * area2;
00968     }
00969     CartVect absnorm( fabs( normal[0] ), fabs( normal[1] ), fabs( normal[2] ) );
00970 
00971     // project to the closest coordinate plane so we only have to do this in 2D
00972 
00973     if( absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2] )
00974     {
00975         area2 = determ3( p[0][1], p[0][2], p[1][1], p[1][2], p[2][1], p[2][2] );
00976         if( fabs( area2 ) < tol )
00977         {
00978             areacoord = CartVect( -std::numeric_limits< double >::min() );  // .set(
00979                                                                             // -std::numeric_limits<double>::min(),
00980                                                                             // -std::numeric_limits<double>::min(),
00981                                                                             // -std::numeric_limits<double>::min() );
00982         }
00983         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
00984         {
00985             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
00986         }
00987         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
00988         {
00989             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
00990         }
00991         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
00992         {
00993             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
00994         }
00995         else
00996         {
00997 
00998             areacoord[0] = determ3( pt_on_plane[1], pt_on_plane[2], p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
00999 
01000             areacoord[1] = determ3( p[0][1], p[0][2], pt_on_plane[1], pt_on_plane[2], p[2][1], p[2][2] ) / area2;
01001 
01002             areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2], pt_on_plane[1], pt_on_plane[2] ) / area2;
01003         }
01004     }
01005     else if( absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2] )
01006     {
01007 
01008         area2 = determ3( p[0][0], p[0][2], p[1][0], p[1][2], p[2][0], p[2][2] );
01009         if( fabs( area2 ) < tol )
01010         {
01011             areacoord = CartVect( -std::numeric_limits< double >::min() );  //.set(
01012                                                                             //-std::numeric_limits<double>::min(),
01013                                                                             //-std::numeric_limits<double>::min(),
01014                                                                             //-std::numeric_limits<double>::min() );
01015         }
01016         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
01017         {
01018             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
01019         }
01020         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
01021         {
01022             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
01023         }
01024         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
01025         {
01026             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
01027         }
01028         else
01029         {
01030 
01031             areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[2], p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
01032 
01033             areacoord[1] = determ3( p[0][0], p[0][2], pt_on_plane[0], pt_on_plane[2], p[2][0], p[2][2] ) / area2;
01034 
01035             areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2], pt_on_plane[0], pt_on_plane[2] ) / area2;
01036         }
01037     }
01038     else
01039     {
01040         /*area2 = determ3(pt0->x(), pt0->y(),
01041          pt1->x(), pt1->y(),
01042          pt2->x(), pt2->y());*/
01043         area2 = determ3( p[0][0], p[0][1], p[1][0], p[1][1], p[2][0], p[2][1] );
01044         if( fabs( area2 ) < tol )
01045         {
01046             areacoord = CartVect( -std::numeric_limits< double >::min() );  //.set(
01047                                                                             //-std::numeric_limits<double>::min(),
01048                                                                             //-std::numeric_limits<double>::min(),
01049                                                                             //-std::numeric_limits<double>::min() );
01050         }
01051         else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
01052         {
01053             areacoord = CartVect( 1., 0., 0. );  //.set( 1.0, 0.0, 0.0 );
01054         }
01055         else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
01056         {
01057             areacoord = CartVect( 0., 1., 0. );  //.set( 0.0, 1.0, 0.0 );
01058         }
01059         else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
01060         {
01061             areacoord = CartVect( 0., 0., 1. );  //.set( 0.0, 0.0, 1.0 );
01062         }
01063         else
01064         {
01065 
01066             areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[1], p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
01067 
01068             areacoord[1] = determ3( p[0][0], p[0][1], pt_on_plane[0], pt_on_plane[1], p[2][0], p[2][1] ) / area2;
01069 
01070             areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1], pt_on_plane[0], pt_on_plane[1] ) / area2;
01071         }
01072     }
01073 }
01074 
01075 ErrorCode SmoothFace::project_to_facets_main( CartVect& this_point,
01076                                               bool trim,
01077                                               bool& outside,
01078                                               CartVect* closest_point_ptr,
01079                                               CartVect* normal_ptr )
01080 {
01081 
01082     // if there are a lot of facets on this surface - use the OBB search first
01083     // to narrow the selection
01084 
01085     _evaluationsCounter++;
01086     double tolerance = 1.e-5;
01087     std::vector< EntityHandle > facets_out;
01088     // we will start with a list of facets anyway, the best among them wins
01089     ErrorCode rval =
01090         _my_geomTopoTool->obb_tree()->closest_to_location( (double*)&this_point, _obb_root, tolerance, facets_out );
01091     if( MB_SUCCESS != rval ) return rval;
01092 
01093     int interpOrder        = 4;
01094     double compareTol      = 1.e-5;
01095     EntityHandle lastFacet = facets_out.front();
01096     rval = project_to_facets( facets_out, lastFacet, interpOrder, compareTol, this_point, trim, outside,
01097                               closest_point_ptr, normal_ptr );
01098 
01099     return rval;
01100 }
01101 ErrorCode SmoothFace::project_to_facets( std::vector< EntityHandle >& facet_list,
01102                                          EntityHandle& lastFacet,
01103                                          int interpOrder,
01104                                          double compareTol,
01105                                          CartVect& this_point,
01106                                          bool,
01107                                          bool& outside,
01108                                          CartVect* closest_point_ptr,
01109                                          CartVect* normal_ptr )
01110 {
01111 
01112     bool outside_facet      = false;
01113     bool best_outside_facet = true;
01114     double mindist          = 1.e20;
01115     CartVect close_point, best_point( mindist, mindist, mindist ), best_areacoord;
01116     EntityHandle best_facet = 0L;  // no best facet found yet
01117     EntityHandle facet;
01118     assert( facet_list.size() > 0 );
01119 
01120     double big_dist = compareTol * 1.0e3;
01121 
01122     // from the list of close facets, determine the closest point
01123     for( size_t i = 0; i < facet_list.size(); i++ )
01124     {
01125         facet = facet_list[i];
01126         CartVect pt_on_plane;
01127         double dist_to_plane;
01128         project_to_facet_plane( facet, this_point, pt_on_plane, dist_to_plane );
01129 
01130         CartVect areacoord;
01131         // CartVect close_point;
01132         facet_area_coordinate( facet, pt_on_plane, areacoord );
01133         if( interpOrder != 0 )
01134         {
01135 
01136             // modify the areacoord - project to the bezier patch- snaps to the
01137             // edge of the patch if necessary
01138 
01139             if( project_to_facet( facet, this_point, areacoord, close_point, outside_facet, compareTol ) != MB_SUCCESS )
01140             {
01141                 return MB_FAILURE;
01142             }
01143             // if (closest_point_ptr)
01144             //*closest_point_ptr = close_point;
01145         }
01146         // keep track of the minimum distance
01147 
01148         double dist = ( close_point - this_point ).length();  // close_point.distance_between(this_point);
01149         if( ( best_outside_facet == outside_facet && dist < mindist ) ||
01150             ( best_outside_facet && !outside_facet && ( dist < big_dist || best_facet == 0L /*!best_facet*/ ) ) )
01151         {
01152             mindist            = dist;
01153             best_point         = close_point;
01154             best_facet         = facet;
01155             best_areacoord     = areacoord;
01156             best_outside_facet = outside_facet;
01157 
01158             if( dist < compareTol )
01159             {
01160                 break;
01161             }
01162             big_dist = 10.0 * mindist;
01163         }
01164         // facet->marked(1);
01165         // used_facet_list.append(facet);
01166     }
01167 
01168     if( normal_ptr )
01169     {
01170         CartVect normal;
01171         if( eval_bezier_patch_normal( best_facet, best_areacoord, normal ) != MB_SUCCESS )
01172         {
01173             return MB_FAILURE;
01174         }
01175         *normal_ptr = normal;
01176     }
01177 
01178     if( closest_point_ptr )
01179     {
01180         *closest_point_ptr = best_point;
01181     }
01182 
01183     outside   = best_outside_facet;
01184     lastFacet = best_facet;
01185 
01186     return MB_SUCCESS;
01187     // end copy
01188 }
01189 
01190 //===========================================================================
01191 // Function Name: project_to_patch
01192 //
01193 // Member Type:  PUBLIC
01194 // Description:  Project a point to a bezier patch. Pass in the areacoord
01195 //              of the point projected to the linear facet.  Function
01196 //              assumes that the point is contained within the patch -
01197 //              if not, it will project to one of its edges.
01198 //===========================================================================
01199 ErrorCode SmoothFace::project_to_patch( EntityHandle facet,   // (IN) the facet where the patch is defined
01200                                         CartVect& ac,         // (IN) area coordinate initial guess (from linear facet)
01201                                         CartVect& pt,         // (IN) point we are projecting to patch
01202                                         CartVect& eval_pt,    // (OUT) The projected point
01203                                         CartVect* eval_norm,  // (OUT) normal at evaluated point
01204                                         bool& outside,        // (OUT) the closest point on patch to pt is on an edge
01205                                         double compare_tol,   // (IN) comparison tolerance
01206                                         int edge_id )         // (IN) only used if this is to be projected to one
01207 //      of the edges.  Otherwise, should be -1
01208 {
01209     ErrorCode status = MB_SUCCESS;
01210 
01211     // see if we are at a vertex
01212 
01213 #define INCR 0.01
01214     const double tol = compare_tol;
01215 
01216     if( is_at_vertex( facet, pt, ac, compare_tol, eval_pt, eval_norm ) )
01217     {
01218         outside = false;
01219         return MB_SUCCESS;
01220     }
01221 
01222     // check if the start ac is inside the patch -if not, then move it there
01223 
01224     int nout          = 0;
01225     const double atol = 0.001;
01226     if( move_ac_inside( ac, atol ) ) nout++;
01227 
01228     int diverge = 0;
01229     int iter    = 0;
01230     CartVect newpt;
01231     eval_bezier_patch( facet, ac, newpt );
01232     CartVect move   = pt - newpt;
01233     double lastdist = move.length();
01234     double bestdist = lastdist;
01235     CartVect bestac = ac;
01236     CartVect bestpt = newpt;
01237     CartVect bestnorm( 0, 0, 0 );
01238 
01239     // If we are already close enough, then return now
01240 
01241     if( lastdist <= tol && !eval_norm && nout == 0 )
01242     {
01243         eval_pt = pt;
01244         outside = false;
01245         return status;
01246     }
01247 
01248     double ratio, mag, umove, vmove, det, distnew, movedist;
01249     CartVect lastpt = newpt;
01250     CartVect lastac = ac;
01251     CartVect norm;
01252     CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
01253     CartVect du, dv, newac;
01254     bool done = false;
01255     while( !done )
01256     {
01257 
01258         // We will be locating the projected point within the u,v,w coordinate
01259         // system of the triangular bezier patch.  Since u+v+w=1, the components
01260         // are not linearly independent.  We will choose only two of the
01261         // coordinates to use and compute the third.
01262 
01263         int system;
01264         if( lastac[0] >= lastac[1] && lastac[0] >= lastac[2] )
01265         {
01266             system = 0;
01267         }
01268         else if( lastac[1] >= lastac[2] )
01269         {
01270             system = 1;
01271         }
01272         else
01273         {
01274             system = 2;
01275         }
01276 
01277         // compute the surface derivatives with respect to each
01278         // of the barycentric coordinates
01279 
01280         if( system == 1 || system == 2 )
01281         {
01282             xac[0] = lastac[0] + INCR;  // xac.x( lastac.x() + INCR );
01283             if( lastac[1] + lastac[2] == 0.0 ) return MB_FAILURE;
01284             ratio  = lastac[2] / ( lastac[1] + lastac[2] );  // ratio = lastac.z() / (lastac.y() + lastac.z());
01285             xac[1] = ( 1.0 - xac[0] ) * ( 1.0 - ratio );     // xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
01286             xac[2] = 1.0 - xac[0] - xac[1];                  // xac.z( 1.0 - xac.x() - xac.y() );
01287             eval_bezier_patch( facet, xac, xpt );
01288             xvec = xpt - lastpt;
01289             xvec /= INCR;
01290         }
01291         if( system == 0 || system == 2 )
01292         {
01293             yac[1] = ( lastac[1] + INCR );      // yac.y( lastac.y() + INCR );
01294             if( lastac[0] + lastac[2] == 0.0 )  // if (lastac.x() + lastac.z() == 0.0)
01295                 return MB_FAILURE;
01296             ratio  = lastac[2] / ( lastac[0] + lastac[2] );   // ratio = lastac.z() / (lastac.x() + lastac.z());
01297             yac[0] = ( ( 1.0 - yac[1] ) * ( 1.0 - ratio ) );  // yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
01298             yac[2] = ( 1.0 - yac[0] - yac[1] );               // yac.z( 1.0 - yac.x() - yac.y() );
01299             eval_bezier_patch( facet, yac, ypt );
01300             yvec = ypt - lastpt;
01301             yvec /= INCR;
01302         }
01303         if( system == 0 || system == 1 )
01304         {
01305             zac[2] = ( lastac[2] + INCR );      // zac.z( lastac.z() + INCR );
01306             if( lastac[0] + lastac[1] == 0.0 )  // if (lastac.x() + lastac.y() == 0.0)
01307                 return MB_FAILURE;
01308             ratio  = lastac[1] / ( lastac[0] + lastac[1] );   // ratio = lastac.y() / (lastac.x() + lastac.y());
01309             zac[0] = ( ( 1.0 - zac[2] ) * ( 1.0 - ratio ) );  // zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
01310             zac[1] = ( 1.0 - zac[0] - zac[2] );               // zac.y( 1.0 - zac.x() - zac.z() );
01311             eval_bezier_patch( facet, zac, zpt );
01312             zvec = zpt - lastpt;
01313             zvec /= INCR;
01314         }
01315 
01316         // compute the surface normal
01317 
01318         switch( system )
01319         {
01320             case 0:
01321                 du = yvec;
01322                 dv = zvec;
01323                 break;
01324             case 1:
01325                 du = zvec;
01326                 dv = xvec;
01327                 break;
01328             case 2:
01329                 du = xvec;
01330                 dv = yvec;
01331                 break;
01332         }
01333         norm = du * dv;
01334         mag  = norm.length();
01335         if( mag < DBL_EPSILON )
01336         {
01337             return MB_FAILURE;
01338             // do something else here (it is likely a flat triangle -
01339             // so try evaluating just an edge of the bezier patch)
01340         }
01341         norm /= mag;
01342         if( iter == 0 ) bestnorm = norm;
01343 
01344         // project the move vector to the tangent plane
01345 
01346         move = ( norm * move ) * norm;
01347 
01348         // compute an equivalent u-v-w vector
01349 
01350         CartVect absnorm( fabs( norm[0] ), fabs( norm[1] ), fabs( norm[2] ) );
01351         if( absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0] )
01352         {
01353             det = du[0] * dv[1] - dv[0] * du[1];
01354             if( fabs( det ) <= DBL_EPSILON )
01355             {
01356                 return MB_FAILURE;  // do something else here
01357             }
01358             umove = ( move[0] * dv[1] - dv[0] * move[1] ) / det;
01359             vmove = ( du[0] * move[1] - move[0] * du[1] ) / det;
01360         }
01361         else if( absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0] )
01362         {
01363             det = du[0] * dv[2] - dv[0] * du[2];
01364             if( fabs( det ) <= DBL_EPSILON )
01365             {
01366                 return MB_FAILURE;
01367             }
01368             umove = ( move[0] * dv[2] - dv[0] * move[2] ) / det;
01369             vmove = ( du[0] * move[2] - move[0] * du[2] ) / det;
01370         }
01371         else
01372         {
01373             det = du[1] * dv[2] - dv[1] * du[2];
01374             if( fabs( det ) <= DBL_EPSILON )
01375             {
01376                 return MB_FAILURE;
01377             }
01378             umove = ( move[1] * dv[2] - dv[1] * move[2] ) / det;
01379             vmove = ( du[1] * move[2] - move[1] * du[2] ) / det;
01380         }
01381 
01382         /* === compute the new u-v coords and evaluate surface at new location */
01383 
01384         switch( system )
01385         {
01386             case 0:
01387                 newac[1] = ( lastac[1] + umove );          // newac.y( lastac.y() + umove );
01388                 newac[2] = ( lastac[2] + vmove );          // newac.z( lastac.z() + vmove );
01389                 newac[0] = ( 1.0 - newac[1] - newac[2] );  // newac.x( 1.0 - newac.y() - newac.z() );
01390                 break;
01391             case 1:
01392                 newac[2] = ( lastac[2] + umove );          // newac.z( lastac.z() + umove );
01393                 newac[0] = ( lastac[0] + vmove );          // newac.x( lastac.x() + vmove );
01394                 newac[1] = ( 1.0 - newac[2] - newac[0] );  // newac.y( 1.0 - newac.z() - newac.x() );
01395                 break;
01396             case 2:
01397                 newac[0] = ( lastac[0] + umove );          // newac.x( lastac.x() + umove );
01398                 newac[1] = ( lastac[1] + vmove );          // newac.y( lastac.y() + vmove );
01399                 newac[2] = ( 1.0 - newac[0] - newac[1] );  // newac.z( 1.0 - newac.x() - newac.y() );
01400                 break;
01401         }
01402 
01403         // Keep it inside the patch
01404 
01405         if( newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol )
01406         {
01407             nout = 0;
01408         }
01409         else
01410         {
01411             if( move_ac_inside( newac, atol ) ) nout++;
01412         }
01413 
01414         // Evaluate at the new location
01415 
01416         if( edge_id != -1 ) ac_at_edge( newac, newac, edge_id );  // move to edge first
01417         eval_bezier_patch( facet, newac, newpt );
01418 
01419         // Check for convergence
01420 
01421         distnew  = ( pt - newpt ).length();  // pt.distance_between(newpt);
01422         move     = newpt - lastpt;
01423         movedist = move.length();
01424         if( movedist < tol || distnew < tol )
01425         {
01426             done = true;
01427             if( distnew < bestdist )
01428             {
01429                 bestdist = distnew;
01430                 bestac   = newac;
01431                 bestpt   = newpt;
01432                 bestnorm = norm;
01433             }
01434         }
01435         else
01436         {
01437 
01438             // don't allow more than 30 iterations
01439 
01440             iter++;
01441             if( iter > 30 )
01442             {
01443                 // if (movedist > tol * 100.0) nout=1;
01444                 done = true;
01445             }
01446 
01447             // Check for divergence - don't allow more than 5 divergent
01448             // iterations
01449 
01450             if( distnew > lastdist )
01451             {
01452                 diverge++;
01453                 if( diverge > 10 )
01454                 {
01455                     done = true;
01456                     // if (movedist > tol * 100.0) nout=1;
01457                 }
01458             }
01459 
01460             // Check if we are continuing to project outside the facet.
01461             // If so, then stop now
01462 
01463             if( nout > 3 )
01464             {
01465                 done = true;
01466             }
01467 
01468             // set up for next iteration
01469 
01470             if( !done )
01471             {
01472                 if( distnew < bestdist )
01473                 {
01474                     bestdist = distnew;
01475                     bestac   = newac;
01476                     bestpt   = newpt;
01477                     bestnorm = norm;
01478                 }
01479                 lastdist = distnew;
01480                 lastpt   = newpt;
01481                 lastac   = newac;
01482                 move     = pt - lastpt;
01483             }
01484         }
01485     }
01486 
01487     eval_pt = bestpt;
01488     if( eval_norm )
01489     {
01490         *eval_norm = bestnorm;
01491     }
01492     outside = ( nout > 0 ) ? true : false;
01493     ac      = bestac;
01494 
01495     return status;
01496 }
01497 
01498 //===========================================================================
01499 // Function Name: ac_at_edge
01500 //
01501 // Member Type:  PRIVATE
01502 // Description:  determine the area coordinate of the facet at the edge
01503 //===========================================================================
01504 void SmoothFace::ac_at_edge( CartVect& fac,  // facet area coordinate
01505                              CartVect& eac,  // edge area coordinate
01506                              int edge_id )   // id of edge
01507 {
01508     double u, v, w;
01509     switch( edge_id )
01510     {
01511         case 0:
01512             u = 0.0;
01513             v = fac[1] / ( fac[1] + fac[2] );  // v = fac.y() / (fac.y() + fac.z());
01514             w = 1.0 - v;
01515             break;
01516         case 1:
01517             u = fac[0] / ( fac[0] + fac[2] );  // u = fac.x() / (fac.x() + fac.z());
01518             v = 0.0;
01519             w = 1.0 - u;
01520             break;
01521         case 2:
01522             u = fac[0] / ( fac[0] + fac[1] );  // u = fac.x() / (fac.x() + fac.y());
01523             v = 1.0 - u;
01524             w = 0.0;
01525             break;
01526         default:
01527             assert( 0 );
01528             u = -1;  // needed to eliminate warnings about used before set
01529             v = -1;  // needed to eliminate warnings about used before set
01530             w = -1;  // needed to eliminate warnings about used before set
01531             break;
01532     }
01533     eac[0] = u;
01534     eac[1] = v;
01535     eac[2] = w;  //= CartVect(u, v, w);
01536 }
01537 
01538 //===========================================================================
01539 // Function Name: project_to_facet
01540 //
01541 // Member Type:  PUBLIC
01542 // Description:  project to a single facet.  Uses the input areacoord as
01543 //              a starting guess.
01544 //===========================================================================
01545 ErrorCode SmoothFace::project_to_facet( EntityHandle facet,
01546                                         CartVect& pt,
01547                                         CartVect& areacoord,
01548                                         CartVect& close_point,
01549                                         bool& outside_facet,
01550                                         double compare_tol )
01551 {
01552     const EntityHandle* conn3 = NULL;
01553     int nnodes                = 0;
01554     _mb->get_connectivity( facet, conn3, nnodes );
01555     //
01556     // double coords[9]; // store the coordinates for the nodes
01557     //_mb->get_coords(conn3, 3, coords);
01558     CartVect p[3];
01559     _mb->get_coords( conn3, 3, (double*)&p[0] );
01560 
01561     int edge_id    = -1;
01562     ErrorCode stat = project_to_patch( facet, areacoord, pt, close_point, NULL, outside_facet, compare_tol, edge_id );
01563     /* }
01564      break;
01565      }*/
01566 
01567     return stat;
01568 }
01569 
01570 //===========================================================================
01571 // Function Name: is_at_vertex
01572 //
01573 // Member Type:  PRIVATE
01574 // Description:  determine if the point is at one of the facet's vertices
01575 //===========================================================================
01576 bool SmoothFace::is_at_vertex( EntityHandle facet,        // (IN) facet we are evaluating
01577                                CartVect& pt,              // (IN) the point
01578                                CartVect& ac,              // (IN) the ac of the point on the facet plane
01579                                double compare_tol,        // (IN) return TRUE of closer than this
01580                                CartVect& eval_pt,         // (OUT) location at vertex if TRUE
01581                                CartVect* eval_norm_ptr )  // (OUT) normal at vertex if TRUE
01582 {
01583     double dist;
01584     CartVect vert_loc;
01585     const double actol = 0.1;
01586     // get coordinates get_coords
01587     const EntityHandle* conn3 = NULL;
01588     int nnodes                = 0;
01589     _mb->get_connectivity( facet, conn3, nnodes );
01590     //
01591     // double coords[9]; // store the coordinates for the nodes
01592     //_mb->get_coords(conn3, 3, coords);
01593     CartVect p[3];
01594     _mb->get_coords( conn3, 3, (double*)&p[0] );
01595     // also get the normals at nodes
01596     CartVect NN[3];
01597     _mb->tag_get_data( _gradientTag, conn3, 3, (double*)&NN[0] );
01598     if( fabs( ac[0] ) < actol && fabs( ac[1] ) < actol )
01599     {
01600         vert_loc = p[2];
01601         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01602         if( dist <= compare_tol )
01603         {
01604             eval_pt = vert_loc;
01605             if( eval_norm_ptr )
01606             {
01607                 *eval_norm_ptr = NN[2];
01608             }
01609             return true;
01610         }
01611     }
01612 
01613     if( fabs( ac[0] ) < actol && fabs( ac[2] ) < actol )
01614     {
01615         vert_loc = p[1];
01616         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01617         if( dist <= compare_tol )
01618         {
01619             eval_pt = vert_loc;
01620             if( eval_norm_ptr )
01621             {
01622                 *eval_norm_ptr = NN[1];  // facet->point(1)->normal( facet );
01623             }
01624             return true;
01625         }
01626     }
01627 
01628     if( fabs( ac[1] ) < actol && fabs( ac[2] ) < actol )
01629     {
01630         vert_loc = p[0];
01631         dist     = ( pt - vert_loc ).length();  // pt.distance_between( vert_loc );
01632         if( dist <= compare_tol )
01633         {
01634             eval_pt = vert_loc;
01635             if( eval_norm_ptr )
01636             {
01637                 *eval_norm_ptr = NN[0];
01638             }
01639             return true;
01640         }
01641     }
01642 
01643     return false;
01644 }
01645 
01646 //===========================================================================
01647 // Function Name: move_ac_inside
01648 //
01649 // Member Type:  PRIVATE
01650 // Description:  find the closest area coordinate to the boundary of the
01651 //              patch if any of its components are < 0
01652 //              Return if the ac was modified.
01653 //===========================================================================
01654 bool SmoothFace::move_ac_inside( CartVect& ac, double tol )
01655 {
01656     int nout = 0;
01657     if( ac[0] < -tol )
01658     {
01659         ac[0] = 0.0;
01660         ac[1] = ac[1] / ( ac[1] + ac[2] );  //( ac.y() / (ac.y() + ac.z()) ;
01661         ac[2] = 1. - ac[1];                 // ac.z( 1.0 - ac.y() );
01662         nout++;
01663     }
01664     if( ac[1] < -tol )
01665     {
01666         ac[1] = 0.;                         // ac.y( 0.0 );
01667         ac[0] = ac[0] / ( ac[0] + ac[2] );  // ac.x( ac.x() / (ac.x() + ac.z()) );
01668         ac[2] = 1. - ac[0];                 // ac.z( 1.0 - ac.x() );
01669         nout++;
01670     }
01671     if( ac[2] < -tol )
01672     {
01673         ac[2] = 0.;                         // ac.z( 0.0 );
01674         ac[0] = ac[0] / ( ac[0] + ac[1] );  // ac.x( ac.x() / (ac.x() + ac.y()) );
01675         ac[1] = 1. - ac[0];                 // ac.y( 1.0 - ac.x() );
01676         nout++;
01677     }
01678     return ( nout > 0 ) ? true : false;
01679 }
01680 
01681 //===========================================================================
01682 // Function Name: hodograph
01683 //
01684 // Member Type:  PUBLIC
01685 // Description:  get the hodograph control points for the facet
01686 // Note:  This is a triangle cubic patch that is defined by the
01687 //       normals of quartic facet control point lattice.  Returned coordinates
01688 //       in Nijk are defined by the following diagram
01689 //
01690 //
01691 //                         *9               index  polar
01692 //                        / \                 0     300    point(0)
01693 //                       /   \                1     210
01694 //                     7*-----*8              2     120
01695 //                     / \   / \              3     030    point(1)
01696 //                    /   \ /   \             4     201
01697 //                  4*----5*-----*6           5     111
01698 //                  / \   / \   / \           6     021
01699 //                 /   \ /   \ /   \          7     102
01700 //                *-----*-----*-----*         8     012
01701 //                0     1     2     3         9     003    point(2)
01702 //
01703 
01704 //===========================================================================
01705 // Function Name: eval_bezier_patch_normal
01706 //
01707 // Member Type:  PRIVATE
01708 // Description:  evaluate the Bezier patch defined at a facet
01709 //===========================================================================
01710 ErrorCode SmoothFace::eval_bezier_patch_normal( EntityHandle facet, CartVect& areacoord, CartVect& normal )
01711 {
01712     // interpolate internal control points
01713 
01714     CartVect gctrl_pts[6];
01715     // facet->get_control_points( gctrl_pts );
01716     ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &facet, 1, &gctrl_pts[0] );
01717     assert( rval == MB_SUCCESS );
01718     if( MB_SUCCESS != rval ) return rval;
01719     // _gradientTag
01720     // get normals at points
01721     const EntityHandle* conn3 = NULL;
01722     int nnodes                = 0;
01723     rval                      = _mb->get_connectivity( facet, conn3, nnodes );
01724     if( MB_SUCCESS != rval ) return rval;
01725 
01726     CartVect NN[3];
01727     rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
01728     assert( rval == MB_SUCCESS );
01729     if( MB_SUCCESS != rval ) return rval;
01730 
01731     if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
01732     {
01733         normal = NN[0];
01734         return MB_SUCCESS;
01735     }
01736     if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
01737     {
01738         normal = NN[1];  // facet->point(1)->normal(facet);
01739         return MB_SUCCESS;
01740     }
01741     if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
01742     {
01743         normal = NN[2];  // facet->point(2)->normal(facet);
01744         return MB_SUCCESS;
01745     }
01746 
01747     // compute the hodograph of the quartic Gregory patch
01748 
01749     CartVect Nijk[10];
01750     // hodograph(facet,areacoord,Nijk);
01751     // start copy from hodograph
01752     // CubitVector gctrl_pts[6];
01753     // facet->get_control_points( gctrl_pts );
01754     CartVect P_facet[3];
01755 
01756     // 2,1,1
01757     /*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
01758      (areacoord.y() * gctrl_pts[3] +
01759      areacoord.z() * gctrl_pts[4]);*/
01760     P_facet[0] =
01761         ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
01762     // 1,2,1
01763     /*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
01764      (areacoord.x() * gctrl_pts[0] +
01765      areacoord.z() * gctrl_pts[5]);*/
01766     P_facet[1] =
01767         ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
01768     // 1,1,2
01769     /*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
01770      (areacoord.x() * gctrl_pts[1] +
01771      areacoord.y() * gctrl_pts[2]);*/
01772     P_facet[2] =
01773         ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
01774 
01775     // corner control points are just the normals at the points
01776 
01777     // 3, 0, 0
01778     Nijk[0] = NN[0];
01779     // 0, 3, 0
01780     Nijk[3] = NN[1];
01781     // 0, 0, 3
01782     Nijk[9] = NN[2];  // facet->point(2)->normal(facet);
01783 
01784     // fill in the boundary control points.  Define as the normal to the local
01785     // triangle formed by the quartic control point lattice
01786 
01787     // store here again the 9 control points on the edges
01788     CartVect CP[9];  // 9 control points on the edges,
01789     rval = _mb->tag_get_data( _facetEdgeCtrlTag, &facet, 1, &CP[0] );
01790     if( MB_SUCCESS != rval ) return rval;
01791     // there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
01792     // CubitFacetEdge *edge;
01793     // edge = facet->edge( 2 );
01794     // CubitVector ctrl_pts[5];
01795     // edge->control_points(facet, ctrl_pts);
01796 
01797     // 2, 1, 0
01798     // Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
01799     Nijk[1] = ( CP[7] - CP[6] ) * ( P_facet[0] - CP[6] );
01800     Nijk[1].normalize();
01801 
01802     // 1, 2, 0
01803     // Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
01804     Nijk[2] = ( CP[8] - CP[7] ) * ( P_facet[1] - CP[7] );
01805     Nijk[2].normalize();
01806 
01807     // edge = facet->edge( 0 );
01808     // edge->control_points(facet, ctrl_pts);
01809 
01810     // 0, 2, 1
01811     // Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
01812     Nijk[6] = ( CP[0] - P_facet[1] ) * ( CP[1] - P_facet[1] );
01813     Nijk[6].normalize();
01814 
01815     // 0, 1, 2
01816     // Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
01817     Nijk[8] = ( CP[1] - P_facet[2] ) * ( CP[2] - P_facet[2] );
01818     Nijk[8].normalize();
01819 
01820     // edge = facet->edge( 1 );
01821     // edge->control_points(facet, ctrl_pts);
01822 
01823     // 1, 0, 2
01824     // Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
01825     Nijk[7] = ( P_facet[2] - CP[4] ) * ( CP[3] - CP[4] );
01826     Nijk[7].normalize();
01827 
01828     // 2, 0, 1
01829     // Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
01830     Nijk[4] = ( P_facet[0] - CP[5] ) * ( CP[4] - CP[5] );
01831     Nijk[4].normalize();
01832 
01833     // 1, 1, 1
01834     Nijk[5] = ( P_facet[1] - P_facet[0] ) * ( P_facet[2] - P_facet[0] );
01835     Nijk[5].normalize();
01836     // end copy from hodograph
01837 
01838     // sum the contribution from each of the control points
01839 
01840     normal = CartVect( 0.0e0, 0.0e0, 0.0e0 );
01841 
01842     // i=3; j=0; k=0;
01843     // double Bsum = 0.0;
01844     double B = mbcube( areacoord[0] );
01845     // Bsum += B;
01846     normal += B * Nijk[0];
01847 
01848     // i=2; j=1; k=0;
01849     B = 3.0 * mbsqr( areacoord[0] ) * areacoord[1];
01850     // Bsum += B;
01851     normal += B * Nijk[1];
01852 
01853     // i=1; j=2; k=0;
01854     B = 3.0 * areacoord[0] * mbsqr( areacoord[1] );
01855     // Bsum += B;
01856     normal += B * Nijk[2];
01857 
01858     // i=0; j=3; k=0;
01859     B = mbcube( areacoord[1] );
01860     // Bsum += B;
01861     normal += B * Nijk[3];
01862 
01863     // i=2; j=0; k=1;
01864     B = 3.0 * mbsqr( areacoord[0] ) * areacoord[2];
01865     // Bsum += B;
01866     normal += B * Nijk[4];
01867 
01868     // i=1; j=1; k=1;
01869     B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
01870     // Bsum += B;
01871     normal += B * Nijk[5];
01872 
01873     // i=0; j=2; k=1;
01874     B = 3.0 * mbsqr( areacoord[1] ) * areacoord[2];
01875     // Bsum += B;
01876     normal += B * Nijk[6];
01877 
01878     // i=1; j=0; k=2;
01879     B = 3.0 * areacoord[0] * mbsqr( areacoord[2] );
01880     // Bsum += B;
01881     normal += B * Nijk[7];
01882 
01883     // i=0; j=1; k=2;
01884     B = 3.0 * areacoord[1] * mbsqr( areacoord[2] );
01885     // Bsum += B;
01886     normal += B * Nijk[8];
01887 
01888     // i=0; j=0; k=3;
01889     B = mbcube( areacoord[2] );
01890     // Bsum += B;
01891     normal += B * Nijk[9];
01892 
01893     // assert(fabs(Bsum - 1.0) < 1e-9);
01894 
01895     normal.normalize();
01896 
01897     return MB_SUCCESS;
01898 }
01899 
01900 ErrorCode SmoothFace::get_normals_for_vertices( const EntityHandle* conn2, CartVect N[2] )
01901 // this method will be called to retrieve the normals needed in the calculation of control edge
01902 // points..
01903 {
01904     // CartVect N[2];
01905     //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
01906     ErrorCode rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
01907     return rval;
01908 }
01909 
01910 ErrorCode SmoothFace::ray_intersection_correct( EntityHandle,       // (IN) the facet where the patch is defined
01911                                                 CartVect& pt,       // (IN) shoot from
01912                                                 CartVect& ray,      // (IN) ray direction
01913                                                 CartVect& eval_pt,  // (INOUT) The intersection point
01914                                                 double& distance,   // (IN OUT) the new distance
01915                                                 bool& outside )
01916 {
01917     // find a point on the smooth surface
01918     CartVect currentPoint = eval_pt;
01919     int numIter           = 0;
01920     double improvement    = 1.e20;
01921     CartVect diff;
01922     while( numIter++ < 5 && improvement > 0.01 )
01923     {
01924         CartVect newPos;
01925 
01926         bool trim = false;  // is it needed?
01927         outside   = true;
01928         CartVect closestPoint;
01929         CartVect normal;
01930 
01931         ErrorCode rval = project_to_facets_main( currentPoint, trim, outside, &newPos, &normal );
01932         if( MB_SUCCESS != rval ) return rval;
01933         assert( rval == MB_SUCCESS );
01934         diff        = newPos - currentPoint;
01935         improvement = diff.length();
01936         // ( pt + t * ray - closest ) % normal = 0;
01937         // intersect tangent plane that goes through closest point with the direction
01938         // t = normal%(closest-pt) / normal%ray;
01939         double dot = normal % ray;  // if it is 0, get out while we can
01940         if( dot < 0.00001 )
01941         {
01942             // bad convergence, get out, do not modify anything
01943             return MB_SUCCESS;
01944         }
01945         double t     = ( ( newPos - pt ) % normal ) / ( dot );
01946         currentPoint = pt + t * ray;
01947     }
01948     eval_pt  = currentPoint;
01949     diff     = currentPoint - pt;
01950     distance = diff.length();
01951     return MB_SUCCESS;
01952 }
01953 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines