MeshKit
1.0
|
00001 #include "meshkit/EdgeMesher.hpp" 00002 #include "meshkit/Matrix.hpp" 00003 #include "meshkit/MKCore.hpp" 00004 #include "meshkit/ModelEnt.hpp" 00005 #include "meshkit/SizingFunction.hpp" 00006 #include "meshkit/RegisterMeshOp.hpp" 00007 #include "moab/ReadUtilIface.hpp" 00008 #include <vector> 00009 #include <math.h> 00010 00011 namespace MeshKit 00012 { 00013 00014 //---------------------------------------------------------------------------// 00015 //Entity Type initilization for edge meshing 00016 moab::EntityType EdgeMesher_tps[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBMAXTYPE}; 00017 const moab::EntityType* EdgeMesher::output_types() 00018 { return EdgeMesher_tps; } 00019 00020 //---------------------------------------------------------------------------// 00021 // Construction Function for Edge Mesher 00022 EdgeMesher::EdgeMesher(MKCore *mk_core, const MEntVector &me_vec) 00023 : MeshScheme(mk_core, me_vec), schemeType(EQUAL), ratio(1.2) 00024 { 00025 00026 } 00027 #if 0 00028 //---------------------------------------------IBERRCHK(g_err, "Trouble get the adjacent geometric nodes on a surface.");------------------------------// 00029 // measure function: compute the distance between the parametric coordinate 00030 // ustart and the parametric coordinate uend 00031 double EdgeMesher::measure(iGeom::EntityHandle ent, double ustart, double uend) const 00032 { 00033 double umin, umax; 00034 00035 //get the minimal and maximal parametrical coordinates of the edge 00036 iGeom::Error err = mk_core()->igeom_instance()->getEntURange(ent, umin, umax); 00037 IBERRCHK(err, "Trouble getting parameter range for edge."); 00038 00039 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to same parameter umax and umin."); 00040 00041 //compute the distance for edges 00042 double measure; 00043 err = mk_core()->igeom_instance()->measure(&ent, 1, &measure); 00044 00045 IBERRCHK(err, "Trouble getting edge measure."); 00046 00047 return measure * (uend - ustart) / (umax - umin); 00048 } 00049 #endif 00050 //---------------------------------------------------------------------------// 00051 // setup function: set up the number of intervals for edge meshing through the 00052 // sizing function 00053 void EdgeMesher::setup_this() 00054 { 00055 //compute the number of intervals for the associated ModelEnts, from the size set on them 00056 //the sizing function they point to, or a default sizing function 00057 for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) 00058 { 00059 ModelEnt *me = mit->first; 00060 00061 //first check to see whether entity is meshed 00062 if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0) 00063 continue; 00064 00065 SizingFunction *sf = mk_core()->sizing_function(me->sizing_function_index()); 00066 if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT && 00067 mk_core()->sizing_function(0)) 00068 sf = mk_core()->sizing_function(0); 00069 00070 if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT) 00071 { 00072 //no sizing set, just assume default #intervals as 4 00073 me->mesh_intervals(4); 00074 me->interval_firmness(DEFAULT); 00075 } 00076 else 00077 { 00078 //check # intervals first, then size, and just choose for now 00079 if (sf->intervals() > 0) 00080 { 00081 if (me->constrain_even() && sf->intervals()%2) 00082 me -> mesh_intervals(sf->intervals()+1); 00083 else 00084 me -> mesh_intervals(sf->intervals()); 00085 me -> interval_firmness(HARD); 00086 } 00087 else if (sf->size()>0) 00088 { 00089 int intervals = me->measure()/sf->size(); 00090 if (!intervals) intervals++; 00091 if (me->constrain_even() && intervals%2) intervals++; 00092 me->mesh_intervals(intervals); 00093 me->interval_firmness(SOFT); 00094 } 00095 else 00096 throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, "Sizing function for edge had neither positive size nor positive intervals."); 00097 } 00098 } 00099 00100 // now call setup_boundary to treat vertices 00101 // Wrong!!!!!!!!! 00102 // setup_boundary(); 00103 /* this is not enough to ensure that vertex mesher will be called before 00104 "this" edge mesher 00105 the case that fell through the cracks was if the end nodes were already setup 00106 then the this_op[0] would not be retrieved, and not inserted "before" the edge mesher MeshOp 00107 */ 00108 int dim=0; 00109 MeshOp * vm = (MeshOp*) mk_core()->construct_meshop(dim); 00110 00111 for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) 00112 { 00113 ModelEnt *me = mit->first; 00114 MEntVector children; 00115 me->get_adjacencies(0, children); 00116 for (unsigned int i=0; i<children.size(); i++) 00117 if (children[i]->is_meshops_list_empty()) 00118 { 00119 vm->add_modelent(children[i]); 00120 } 00121 } 00122 // in any case, make sure that the vertex mesher is inserted before this edge mesher 00123 mk_core()->insert_node(vm, this, mk_core()->root_node()); 00124 00125 00126 } 00127 00128 //---------------------------------------------------------------------------// 00129 // execute function: Generate the mesh for edges 00130 void EdgeMesher::execute_this() 00131 { 00132 std::vector<double> coords; 00133 std::vector<moab::EntityHandle> nodes; 00134 00135 for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) 00136 { 00137 ModelEnt *me = mit -> first; 00138 if (me->get_meshed_state() >= COMPLETE_MESH) 00139 continue; 00140 //resize the coords based on the interval setting 00141 int num_edges = me->mesh_intervals(); 00142 coords.resize(3*(num_edges+1)); 00143 nodes.clear(); 00144 nodes.reserve(num_edges + 1); 00145 00146 //get bounding mesh entities, use 1st 2 entries of nodes list temporarily 00147 //pick up the boundary end nodes 00148 me->boundary(0, nodes); 00149 00150 bool periodic = (nodes.size() == 1); 00151 00152 //get coords in list, then move one tuple to the last position 00153 moab::ErrorCode rval = mk_core()->moab_instance()->get_coords(&nodes[0], nodes.size(), &coords[0]); 00154 MBERRCHK(rval, mk_core()->moab_instance()); 00155 00156 //move the second node to the endmost position in the node list 00157 // if periodic, the last node coordinates are also at index 0 in coords array 00158 // there is only one node, coords[3+i] are not even initialized 00159 int index2 = (periodic) ? 0 : 3; 00160 for (int i = 0; i < 3; i++) 00161 coords[3*num_edges+i] = coords[index2+i]; 00162 00163 EdgeSchemeType scheme = schemeType; 00164 SizingFunction *sz = mk_core()->sizing_function(me->sizing_function_index()); 00165 if (sz->variable()) 00166 scheme = VARIABLE; 00167 00168 //choose the scheme for edge mesher 00169 switch(scheme) 00170 { 00171 case EQUAL://equal meshing for edges 00172 EqualMeshing(me, num_edges, coords); 00173 break; 00174 case BIAS://bias meshing for edges 00175 BiasMeshing(me, num_edges, coords); 00176 break; 00177 case DUAL://dual bias meshing for edges 00178 DualBiasMeshing(me, num_edges, coords); 00179 break; 00180 case CURVATURE://curvature-based meshing for edges 00181 CurvatureMeshing(me, num_edges, coords); 00182 break; 00183 case VARIABLE: // use a var size from sizing function 00184 VariableMeshing(me, num_edges, coords); 00185 break; 00186 case EQUIGNOMONIC: // used to generate HOMME type meshes on a sphere 00187 EquiAngleGnomonic(me, num_edges, coords); 00188 break; 00189 default: 00190 break; 00191 } 00192 //the variable nodes should be resized, node size may be changed in the different scheme for edge meshing 00193 me->mesh_intervals(num_edges); 00194 nodes.resize(num_edges+1); 00195 00196 //move the other nodes to the end of nodes' list 00197 if (periodic) nodes[num_edges] = nodes[0]; 00198 else nodes[num_edges] = nodes[1]; 00199 00200 //create the vertices' entities on the edge 00201 if (num_edges > 1) { 00202 rval = mk_core()->moab_instance()->create_vertices(&coords[3], num_edges - 1, mit -> second); 00203 MBERRCHK(rval, mk_core()->moab_instance()); 00204 } 00205 00206 //distribute the nodes into vector 00207 int j = 1; 00208 for (moab::Range::iterator rit = mit -> second.begin(); rit != mit -> second.end(); rit++) 00209 nodes[j++] = *rit; 00210 00211 //get the query interface, which we will use to create the edges directly 00212 moab::ReadUtilIface *iface; 00213 rval = mk_core() -> moab_instance() -> query_interface(iface); 00214 MBERRCHK(rval, mk_core()->moab_instance()); 00215 00216 //create the edges, get a direct ptr to connectivity 00217 moab::EntityHandle starth, *connect, *tmp_connect; 00218 rval = iface -> get_element_connect(num_edges, 2, moab::MBEDGE, 1, starth, connect); 00219 MBERRCHK(rval, mk_core()->moab_instance()); 00220 00221 //add edges to range for the MESelection 00222 mit -> second.insert(starth, starth + num_edges - 1); 00223 00224 //now set the connectivity array from the nodes 00225 tmp_connect = &nodes[0]; 00226 for (int i = 0; i < num_edges; i++) 00227 { 00228 connect[0] = tmp_connect[0]; 00229 connect[1] = tmp_connect[1]; 00230 00231 //increment connectivity ptr by 2 to go to next edge 00232 connect += 2; 00233 00234 //increment tmp_connect by 1, to go to next node 00235 tmp_connect++; 00236 } 00237 00238 // ok, we are done, commit to ME 00239 me->commit_mesh(mit->second, COMPLETE_MESH); 00240 00241 00242 } 00243 00244 00245 } 00246 00247 //---------------------------------------------------------------------------// 00248 // Deconstruction function 00249 EdgeMesher::~EdgeMesher() 00250 { 00251 00252 } 00253 00254 //---------------------------------------------------------------------------// 00255 // Create the mesh for edges with the equal distances 00256 void EdgeMesher::EqualMeshing(ModelEnt *ent, int num_edges, std::vector<double> &coords) 00257 { 00258 double umin, umax, measure; 00259 (void) measure; 00260 //get the u range for the edge 00261 iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax); 00262 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00263 00264 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00265 00266 //get the arc length 00267 measure = ent -> measure(); 00268 00269 double u, du; 00270 if (!num_edges) throw Error(MK_BAD_INPUT, "Trying to mesh edge with zero edges."); 00271 du = (umax - umin)/(double)num_edges; 00272 00273 u = umin; 00274 //distribute the nodes with equal distances 00275 for (int i = 1; i < num_edges; i++) 00276 { 00277 u = umin + i*du; 00278 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]); 00279 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00280 } 00281 00282 } 00283 00284 //---------------------------------------------------------------------------// 00285 // Create the mesh for edges based on the curvatures 00286 void EdgeMesher::CurvatureMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords) 00287 { 00288 double umin, umax, measure; 00289 (void) measure; 00290 //store the initial edge size, the edge size may be changed during meshing 00291 int initial_num_edges = num_edges; 00292 00293 //get the u range for the edge 00294 iGeom::Error gerr = ent->igeom_instance() ->getEntURange(ent->geom_handle(), umin, umax); 00295 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00296 00297 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00298 00299 //get the arc length 00300 measure = ent -> measure(); 00301 00302 int index = 0; 00303 double u, du, uMid; 00304 du = (umax - umin)/(double)num_edges; 00305 00306 std::vector<double> NodeCoordinates; 00307 std::vector<double> TempNode; 00308 std::vector<double> URecord; //record the value of U 00309 00310 Point3D pts0, pts1, ptsMid; 00311 double tmp[3]; 00312 00313 NodeCoordinates.resize(3*(num_edges + 1)); 00314 00315 TempNode.resize(3*1); 00316 URecord.resize(1); 00317 00318 gerr = ent -> igeom_instance() -> getEntUtoXYZ(ent->geom_handle(), umin, TempNode[0], TempNode[1], TempNode[2]); 00319 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge"); 00320 00321 NodeCoordinates[3*0] = TempNode[0]; 00322 NodeCoordinates[3*0+1] = TempNode[1]; 00323 NodeCoordinates[3*0+2] = TempNode[2]; 00324 00325 URecord[0] = umin; 00326 00327 u = umin; 00328 //distribute the mesh nodes on the edge based on the curvature 00329 for (int i = 1; i < num_edges; i++) 00330 { 00331 //first distribute the nodes evenly 00332 u = umin + i*du; 00333 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, NodeCoordinates[3*i], NodeCoordinates[3*i+1], NodeCoordinates[3*i+2]); 00334 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00335 00336 //store the two adjacent mesh nodes on the edge 00337 pts0.px = NodeCoordinates[3*(i-1)]; 00338 pts0.py = NodeCoordinates[3*(i-1)+1]; 00339 pts0.pz = NodeCoordinates[3*(i-1)+2]; 00340 00341 pts1.px = NodeCoordinates[3*i]; 00342 pts1.py = NodeCoordinates[3*i+1]; 00343 pts1.pz = NodeCoordinates[3*i+2]; 00344 00345 //get the coordinates for mid point between two adjacent mesh nodes on the edge 00346 uMid = (u-du+u)/2; 00347 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uMid, tmp[0], tmp[1], tmp[2]); 00348 ptsMid.px = tmp[0]; 00349 ptsMid.py = tmp[1]; 00350 ptsMid.pz = tmp[2]; 00351 00352 //calculate the error and check whether it requires further refinement based on the curvature 00353 if (!ErrorCalculate(ent, pts0, pts1, ptsMid)) 00354 { 00355 DivideIntoMore(ent, pts0, ptsMid, pts1, u-du, u, uMid, index, TempNode, URecord); 00356 } 00357 00358 // add the other end node to the array, get the next two adjacent mesh nodes 00359 index++; 00360 TempNode.resize(3*(index + 1)); 00361 URecord.resize(index + 1); 00362 TempNode[3*index] = pts1.px; 00363 TempNode[3*index + 1] = pts1.py; 00364 TempNode[3*index + 2] = pts1.pz; 00365 00366 URecord[index] = u; 00367 } 00368 00369 //sorting the parametrical coordinate data based on the value of u 00370 assert(TempNode.size()== (3*URecord.size()) ); 00371 00372 QuickSorting(TempNode, URecord, URecord.size()); 00373 num_edges = URecord.size() - 1; 00374 00375 //resize the variable coords 00376 coords.resize(3*(num_edges+1)); 00377 00378 //move the other end node to the endmost of the list 00379 for (int i = 0; i < 3; i++) 00380 coords[3*num_edges+i] = coords[3*initial_num_edges+i]; 00381 00382 //return the mesh nodes 00383 for (int i = 1; i < num_edges; i++) 00384 { 00385 coords[3*i] = TempNode[3*i]; 00386 coords[3*i+1] = TempNode[3*i+1]; 00387 coords[3*i+2] = TempNode[3*i+2]; 00388 } 00389 00390 } 00391 00392 //---------------------------------------------------------------------------// 00393 // Create the mesh for edges with dual bias distances 00394 void EdgeMesher::DualBiasMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords) 00395 { 00396 double umin, umax, measure; 00397 00398 //get the u range for the edge 00399 iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax); 00400 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00401 00402 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00403 00404 //get the arc length 00405 measure = ent -> measure(); 00406 00407 double u, L0, dist, u0, u1; 00408 00409 //if the node # is odd, node # will increase by 1 00410 if ((num_edges%2)!=0) 00411 { 00412 num_edges++; 00413 coords.resize(3*(num_edges+1)); 00414 //move the other end node's position because the variable coords has been resized. 00415 for (int k = 0; k < 3; k++) 00416 coords[3*num_edges + k] = coords[3*num_edges + k - 3]; 00417 } 00418 00419 //default bias ratio is 1.2 00420 double q = ratio;// 00421 00422 //get the distance between the first two nodes 00423 L0 = 0.5 * measure * (1-q) / (1 - pow(q, num_edges/2)); 00424 00425 00426 u0 = umin; 00427 u1 = umax; 00428 //distribute the mesh nodes on the edge with dual-bias distances 00429 for (int i = 1; i < (num_edges/2 + 1); i++) 00430 { 00431 //distribute the mesh nodes on the one side 00432 dist = L0*pow(q, i-1); 00433 u = u0 + (umax - umin) * dist/measure; 00434 u = getUCoord(ent, u0, dist, u, umin, umax); 00435 u0 = u; 00436 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]); 00437 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00438 00439 //distribute the mesh nodes on the other side 00440 if (i < num_edges/2) 00441 { 00442 u = u1 - (umax-umin) * dist / measure; 00443 u = getUCoord(ent, u1, dist, u, umin, umax); 00444 u1 = u; 00445 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*(num_edges-i)], coords[3*(num_edges-i)+1], coords[3*(num_edges-i)+2]); 00446 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge"); 00447 } 00448 00449 } 00450 00451 } 00452 00453 //---------------------------------------------------------------------------// 00454 // Create the mesh for edges with bias distances 00455 void EdgeMesher::BiasMeshing(ModelEnt *ent, int num_edges, std::vector<double> &coords) 00456 { 00457 double umin, umax, measure; 00458 00459 //get the u range for the edge 00460 iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax); 00461 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00462 00463 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00464 00465 //get the arc length 00466 measure = ent -> measure(); 00467 00468 double u, L0, dist = 0, u0; 00469 00470 // the default bias ratio 1.2 00471 double q = ratio; 00472 L0 = measure * (1-q) / (1 - pow(q, num_edges)); 00473 00474 00475 u0 = umin; 00476 //distribute the mesh nodes on the edge with bias distances 00477 for (int i = 1; i < num_edges; i++) 00478 { 00479 dist = L0*pow(q, i-1); 00480 u = u0 + (umax - umin)*dist/measure; 00481 u = getUCoord(ent, u0, dist, u, umin, umax); 00482 u0 = u; 00483 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]); 00484 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00485 } 00486 } 00487 00488 //---------------------------------------------------------------------------// 00489 // Mesh Refinement function based on the curvature: if the error is too big, refine the mesh on the edge 00490 void EdgeMesher::DivideIntoMore(ModelEnt *ent, Point3D p0, Point3D pMid, Point3D p1, double u0, double u1, double uMid, int &index, vector<double> &nodes, vector<double> &URecord) 00491 { 00492 //this is a recursive process, the process continues until the error is smaller than what is required. 00493 //first get two adjacent mesh nodes on the edge and coordinates for mid point between two adjacent nodes. 00494 //then check the left side and right side whether the error is too big nor not 00495 double uu0, uu1, uumid, tmp[3]; 00496 Point3D pts0, pts1, ptsMid; 00497 00498 index++; 00499 nodes.resize(3*(index+1)); 00500 URecord.resize(index+1); 00501 nodes[3*index] = pMid.px; 00502 nodes[3*index+1] = pMid.py; 00503 nodes[3*index+2] = pMid.pz; 00504 URecord[index] = uMid; 00505 00506 //check the left side 00507 uu0=u0; 00508 uu1=uMid; 00509 uumid=(uu0+uu1)/2; 00510 pts0=p0; 00511 pts1=pMid; 00512 00513 00514 iGeom::Error gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uumid, tmp[0], tmp[1], tmp[2]); 00515 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00516 ptsMid.px = tmp[0]; 00517 ptsMid.py = tmp[1]; 00518 ptsMid.pz = tmp[2]; 00519 00520 //check the error 00521 if(!ErrorCalculate(ent, pts0, pts1, ptsMid)) 00522 { 00523 //further refinement 00524 DivideIntoMore(ent, pts0, ptsMid, pts1, uu0, uu1, uumid, index, nodes, URecord); 00525 } 00526 00527 //check the right side 00528 uu0 = uMid; 00529 uu1=u1; 00530 uumid=(uu0+uu1)/2; 00531 pts0=pMid; 00532 pts1=p1; 00533 //get the coorindinates for mid point 00534 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uumid, tmp[0], tmp[1], tmp[2]); 00535 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00536 ptsMid.px = tmp[0]; 00537 ptsMid.py = tmp[1]; 00538 ptsMid.pz = tmp[2]; 00539 00540 //check the error 00541 if(!ErrorCalculate(ent, pts0, pts1, ptsMid)) 00542 { 00543 //further refinement 00544 DivideIntoMore(ent, pts0, ptsMid, pts1, uu0, uu1, uumid, index, nodes, URecord); 00545 } 00546 00547 } 00548 00549 //create the mesh for edges based on variable size from SizingFunction (var) 00550 void EdgeMesher::VariableMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords) 00551 { 00552 double umin, umax, measure; 00553 (void) measure; 00554 // because of that, keep track of the first node position and last node position 00555 // first node position does not change, but the last node position do change 00556 // coords will contain all nodes, including umax in Urecord! 00557 00558 SizingFunction *sf = mk_core()->sizing_function(ent->sizing_function_index()); 00559 //get the u range for the edge 00560 iGeom::EntityHandle edge = ent->geom_handle(); 00561 iGeom::Error gerr = ent->igeom_instance() ->getEntURange(edge, umin, umax); 00562 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00563 00564 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00565 00566 //get the arc length 00567 measure = ent -> measure(); 00568 00569 // start advancing for each edge mesh, from the first point position 00570 double currentPar = umin; 00571 double currentPosition[3]; 00572 gerr = ent->igeom_instance() ->getEntUtoXYZ(edge, umin, currentPosition[0], 00573 currentPosition[1], currentPosition[2] ); 00574 00575 double endPoint[3]; 00576 gerr = ent->igeom_instance() ->getEntUtoXYZ(edge, umax, endPoint[0], 00577 endPoint[1], endPoint[2] ); 00578 Vector<3> endpt(endPoint); 00579 00580 double targetSize = sf->size(currentPosition); 00581 double startSize = targetSize; 00582 00583 double endSize = sf->size(endPoint); 00584 // advance u such that the next point is at "currentSize" distance 00585 // or close to it 00586 // try first with a u that is coming from the (umax-umin)/number of edges 00587 double deltaU = (umax-umin)/num_edges; 00588 //coords.clear(); we do not want to clear, as the first node is still fine 00589 std::vector<double> URecord; //record the values for u; we may have to adjust all 00590 // of them accordingly, and even add one more if we have some evenify problems. 00591 // keep in mind that size is just a suggestion, number of intervals is more important for 00592 // paver mesher 00593 Vector<3> pt(currentPosition); 00594 00595 //bool notDone = true; 00596 double prevU = umin; 00597 while (currentPar + 1.1*deltaU < umax) 00598 { 00599 // do some binary search; actually, better, do Newton-Raphson, which should converge 00600 // faster 00601 // 00602 prevU = currentPar; 00603 currentPar += deltaU; 00604 // adjust current par, such that 00605 double point[3]; 00606 gerr=ent->igeom_instance()->getEntUtoXYZ(edge, currentPar, point[0], point[1], point[2] ); 00607 IBERRCHK(gerr, "Trouble getting position at parameter u."); 00608 Vector<3> ptCandidate(point); 00609 double compSize = length(ptCandidate-pt); 00610 int nit = 0; 00611 00612 while ( (fabs(1.-compSize/targetSize)> 0.02 ) && (nit < 5))// 2% of the target size 00613 { 00614 // do Newton iteration 00615 double tangent[3]; 00616 gerr=ent->igeom_instance() ->getEntTgntU(edge, currentPar, tangent[0], tangent[1], tangent[2] ); 00617 IBERRCHK(gerr, "Trouble getting tangent at parameter u."); 00618 Vector<3> tang(tangent); 00619 double dldu = 1./compSize * ((ptCandidate-pt )%tang); 00620 nit++;// increase iteration count 00621 if (dldu!=0.) 00622 { 00623 double deu= (targetSize-compSize)/dldu; 00624 currentPar+=deu; 00625 if (prevU>currentPar) 00626 { 00627 break; // something is wrong 00628 } 00629 if (umax < currentPar) 00630 { 00631 currentPar = umax; 00632 break; 00633 } 00634 ent->igeom_instance()->getEntUtoXYZ(edge, currentPar, point[0], point[1], point[2]); 00635 Vector<3> newPt(point); 00636 compSize = length(newPt-pt); 00637 ptCandidate = newPt; 00638 } 00639 00640 } 00641 // we have found an acceptable point/param 00642 URecord.push_back(currentPar); 00643 deltaU = currentPar-prevU;// should be greater than 0 00644 pt = ptCandidate; 00645 targetSize = sf->size(pt.data());// a new target size, at the current point 00646 00647 00648 00649 } 00650 // when we are out of here, we need to adjust the URecords, to be more nicely 00651 // distributed; also, look at the evenify again 00652 int sizeU = (int)URecord.size(); 00653 if ((sizeU%2==0) && ent->constrain_even() ) 00654 { 00655 // add one more 00656 if (sizeU==0) 00657 { 00658 // just add one in the middle, and call it done 00659 URecord.push_back( (umin+umax)/2); 00660 } 00661 else 00662 { 00663 //at least 2 (maybe 4) 00664 double lastDelta = URecord[sizeU-1]-URecord[sizeU-2]; 00665 URecord.push_back(URecord[sizeU-1]+lastDelta ); 00666 } 00667 } 00668 // now, we have to redistribute, such as the last 2 deltas are about the same 00669 // so, we should have after a little work, 00670 // umin, umin+c*(URecord[0]-umin), ... umin+c*(URecord[size-1]-umin), umax 00671 // what we have now is 00672 // umin, URecord[0], ... ,URecord[size-1], and umax could be even outside or inside 00673 // keep the last sizes equal 00674 // umin+c(UR[size-2]-umin) + umax = 2*( umin+c*(UR[size-1]-umin)) 00675 // c ( 2*(UR[size-1]-umin) -(UR[size-2]-umin) ) = umax - umin 00676 // c ( 2*UR[size-1] - UR[size-2] - umin ) = umax - umin 00677 // c = (umax-umin)/( 2*UR[size-1] - UR[size-2] - umin) 00678 sizeU = (int)URecord.size();// it may be bigger by one than the last time 00679 if (sizeU == 0) 00680 { 00681 // nothing to do, only one edge to generate 00682 } 00683 else if (sizeU == 1) 00684 { 00685 // put it according to the sizes at ends, and assume a nice variation for u 00686 // (u-umin) / (umax-u) = startSize / endSize 00687 // (u-umin)*endSize = (umax-u) * startSize 00688 // u(endSize+startSize)=(umax*startSize+umin*endSize) 00689 URecord[0] = (umax*startSize+umin*endSize)/(startSize+endSize); 00690 00691 } 00692 else // sizeU>=2, so we can spread the param a little more, assuming nice 00693 // uniform mapping 00694 { 00695 double c = (umax-umin)/( 2*URecord[sizeU-1] - URecord[sizeU-2] - umin); 00696 for (int i=0; i<sizeU; i++) 00697 URecord[i] = umin + c*(URecord[i] -umin);// some spreading out 00698 } 00699 // now, we can finally get the points for each U, U's should be spread out nicely 00700 URecord.push_back(umax); // just add the last u, for the end point 00701 // 00702 sizeU = (int) URecord.size(); // new size, after pushing the last u, umax 00703 num_edges = sizeU;// this is the new number of edges; the last one will be the end point 00704 // of the edge, corresponding to umax 00705 coords.resize(3*sizeU+3); 00706 // we already know that at i=0 is the first node, start vertex of edge 00707 // the rest will be computed from u 00708 // even the last one, which is an overkill 00709 for (int i = 1; i <= num_edges; i++) 00710 { 00711 double u = URecord[i-1]; 00712 gerr = ent->igeom_instance()->getEntUtoXYZ(edge, u, coords[3*i], coords[3*i+1], coords[3*i+2]); 00713 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00714 } 00715 return; 00716 00717 } 00718 // number of edges is input here 00719 // equal angles are formed at the center of the sphere/cube mesh 00720 // it is close to bias meshing, but not quite 00721 void EdgeMesher::EquiAngleGnomonic(ModelEnt *ent, int num_edges, std::vector<double> &coords) 00722 { 00723 const double MY_PI=3.14159265; 00724 double deltaAngle=MY_PI/num_edges/2; 00725 // double length=ent->measure();// this is an edge 00726 double umin, umax; 00727 00728 //get the u range for the edge 00729 iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax); 00730 IBERRCHK(gerr, "Trouble get parameter range for edge."); 00731 00732 if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin."); 00733 00734 // consider that the parametrization is very linear 00735 // most of the time u will be from 0 to length of edge, for a cube 00736 double deltau = umax - umin; 00737 00738 double u = umin;// u will get different values, 00739 // start at u 00740 for (int i = 1; i < num_edges; i++) 00741 { 00742 double betak=i*deltaAngle; 00743 double alfak = MY_PI/4-betak; 00744 double tang_alfak = tan(alfak); 00745 u = umin+deltau/2*(1-tang_alfak); 00746 00747 gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]); 00748 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00749 } 00750 return; 00751 } 00752 //---------------------------------------------------------------------------// 00753 // Rapid sorting the mesh nodes on the edge based on the parametric coordinates. This is a recursive 00754 // process 00755 void EdgeMesher::RapidSorting(vector<double> &nodes, vector<double> &URecord, int left, int right) 00756 { 00757 int i, j; 00758 double middle, iTemp; 00759 Point3D TempData; 00760 00761 middle=URecord[(left+right)/2]; 00762 i=left; 00763 j=right; 00764 00765 do 00766 { 00767 //search the values which are greater than the middle value from the left side 00768 while((URecord[i] < middle)&&(i<right)) 00769 { 00770 i++; 00771 } 00772 //search the values which are greater than the middle value from the right side 00773 while((URecord[j] > middle)&&(j > left)) 00774 { 00775 j--; 00776 } 00777 if (i<=j)//find a pair of values 00778 { 00779 iTemp = URecord[i]; 00780 URecord[i] = URecord[j]; 00781 URecord[j]=iTemp; 00782 00783 00784 TempData.px = nodes[3*i]; 00785 TempData.py = nodes[3*i+1]; 00786 TempData.pz = nodes[3*i+2]; 00787 00788 nodes[3*i] = nodes[3*j]; 00789 nodes[3*i+1] = nodes[3*j+1]; 00790 nodes[3*i+2] = nodes[3*j+2]; 00791 nodes[3*j] = TempData.px; 00792 nodes[3*j+1] = TempData.py; 00793 nodes[3*j+2] = TempData.pz; 00794 00795 i++; 00796 j--; 00797 } 00798 }while(i<=j); 00799 if (left < j) 00800 RapidSorting(nodes, URecord, left, j); 00801 if (right > i) 00802 RapidSorting(nodes, URecord, i, right); 00803 } 00804 00805 //---------------------------------------------------------------------------// 00806 // Quick Sorting: this function comes together with the function RapidSorting 00807 void EdgeMesher::QuickSorting(vector<double> &nodes, vector<double> &URecord, int count) 00808 { 00809 RapidSorting(nodes, URecord, 0, count-1); 00810 } 00811 00812 00813 //---------------------------------------------------------------------------// 00814 // return the x, y, z coordinates from the parametric coordinates 00815 EdgeMesher::Point3D EdgeMesher::getXYZCoords(ModelEnt *ent, double u) const 00816 { 00817 Point3D pts3D; 00818 double xyz[3]; 00819 00820 00821 //get the coordinates in the physical space 00822 iGeom::Error gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, xyz[0], xyz[1], xyz[2]); 00823 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00824 00825 pts3D.px = xyz[0]; 00826 pts3D.py = xyz[1]; 00827 pts3D.pz = xyz[2]; 00828 return pts3D; 00829 } 00830 00831 //---------------------------------------------------------------------------// 00832 // get the parametric coordinate based on first parametric coordinate ustart and distance in the physical space 00833 double EdgeMesher::getUCoord(ModelEnt *ent, double ustart, double dist, double uguess, double umin, double umax) const 00834 { 00835 00836 Point3D p0 = getXYZCoords(ent, ustart); 00837 Point3D p1 = getXYZCoords(ent, uguess); 00838 00839 00840 double dx, dy, dz, dl, u=uguess; 00841 double tol = 1.0E-7; 00842 int test=0; 00843 00844 int ntrials=0; 00845 while(1) 00846 { 00847 dx = p1.px - p0.px; 00848 dy = p1.py - p0.py; 00849 dz = p1.pz - p0.pz; 00850 dl = sqrt(dx * dx + dy * dy + dz * dz); 00851 if ( fabs(dl-dist) < tol) break; 00852 00853 u = ustart + (u - ustart) * (dist/dl); 00854 if (u > umax) 00855 { 00856 u=umax; 00857 test++; 00858 if (test>10) break; 00859 } 00860 if (u < umin) 00861 { 00862 u=umin; 00863 test++; 00864 if (test>10) break; 00865 } 00866 p1 = getXYZCoords(ent, u); 00867 00868 00869 if (ntrials++ == 100000) 00870 { 00871 cout << " Warning: Searching for U failed " << endl; 00872 } 00873 } 00874 uguess = u; 00875 return uguess; 00876 } 00877 00878 //---------------------------------------------------------------------------// 00879 // calculate the error: the distance between the mid point of two adjacent 00880 // mesh nodes (on the mesh line segments) and mid point on the edge 00881 bool EdgeMesher::ErrorCalculate(ModelEnt *ent, Point3D p0, Point3D p1, Point3D pMid) 00882 { 00883 double lengtha, lengthb, lengthc; 00884 double deltax, deltay, deltaz; 00885 double angle, error, tol=1.0E-3, H; 00886 double cvtr_ijk[3], curvature; 00887 bool result; 00888 00889 //calculate the distance between the first mesh node and mid point on the edge 00890 deltax = pMid.px-p0.px; 00891 deltay = pMid.py-p0.py; 00892 deltaz = pMid.pz-p0.pz; 00893 lengtha = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz); 00894 00895 //calculate the distance between two adjacent mesh nodes 00896 deltax = p1.px-p0.px; 00897 deltay = p1.py-p0.py; 00898 deltaz = p1.pz-p0.pz; 00899 lengthb = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz); 00900 00901 //calculate the distance between the second mesh node and mid point on the edge 00902 deltax = pMid.px-p1.px; 00903 deltay = pMid.py-p1.py; 00904 deltaz = pMid.pz-p1.pz; 00905 lengthc = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz); 00906 00907 //calculate the angle 00908 angle = acos((lengtha*lengtha + lengthb*lengthb - lengthc*lengthc)/(2*lengtha*lengthb)); 00909 H = fabs(lengtha*sin(angle)); 00910 00911 //calculate the curvature 00912 iGeom::Error gerr = ent->igeom_instance()->getEgCvtrXYZ(ent->geom_handle(), pMid.px, pMid.py, pMid.pz, cvtr_ijk[0], cvtr_ijk[1], cvtr_ijk[2]); 00913 IBERRCHK(gerr, "Trouble getting U from XYZ along the edge."); 00914 curvature = sqrt(cvtr_ijk[0]*cvtr_ijk[0]+cvtr_ijk[1]*cvtr_ijk[1]+cvtr_ijk[2]*cvtr_ijk[2]); 00915 error= H*curvature; 00916 00917 00918 if (error > tol) 00919 result = false; 00920 else 00921 result = true; 00922 return result; 00923 } 00924 00925 } 00926