MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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