moab
|
00001 #include "moab/DrawDual.hpp" 00002 #include "SheetDiagramPopup.h" 00003 #include "moab/MeshTopoUtil.hpp" 00004 #include "MBTagConventions.hpp" 00005 #include "moab/CartVect.hpp" 00006 #include "vtkMOABUtils.h" 00007 #include "vtkPolyData.h" 00008 #include "vtkPolyDataMapper2D.h" 00009 #include "vtkPolyDataMapper.h" 00010 #include "vtkDataSetMapper.h" 00011 #include "vtkMapper2D.h" 00012 #include "vtkMapper.h" 00013 #include "vtkLabeledDataMapper.h" 00014 #include "vtkActor2D.h" 00015 #include "vtkTextActor.h" 00016 #include "vtkIdFilter.h" 00017 #include "vtkTextProperty.h" 00018 #include "vtkPoints.h" 00019 #include "vtkPointData.h" 00020 #include "vtkRenderer.h" 00021 #include "vtkRenderWindow.h" 00022 #include "vtkRendererCollection.h" 00023 #include "vtkActorCollection.h" 00024 #include "vtkCellArray.h" 00025 #include "vtkCellData.h" 00026 #include "vtkCellPicker.h" 00027 #include "vtkCallbackCommand.h" 00028 #include "vtkCamera.h" 00029 #include "vtkCoordinate.h" 00030 #include "vtkProperty2D.h" 00031 #include "vtkProperty.h" 00032 #include "vtkExtractEdges.h" 00033 #include "vtkExtractCells.h" 00034 #include "vtkTubeFilter.h" 00035 #include "vtkGeometryFilter.h" 00036 #include "vtkInteractorStyle.h" 00037 #include "vtkLookupTable.h" 00038 #include "QVTKWidget.h" 00039 #include "qlineedit.h" 00040 #include "vtkFloatArray.h" 00041 #include "vtkMaskFields.h" 00042 #include "vtkWindowToImageFilter.h" 00043 #include "vtkPNGWriter.h" 00044 #include "assert.h" 00045 #include <sstream> 00046 00047 extern "C" 00048 { 00049 #include "gvc.h" 00050 // extern GVC_t *gvContext(); 00051 } 00052 00053 const bool my_debug = false; 00054 00055 const int SHEET_WINDOW_SIZE = 500; 00056 00057 //const int RAD_PTS = 400; 00058 const int RAD_PTS = 10; 00059 const int CENT_X = 0; 00060 const int CENT_Y = 0; 00061 00062 using namespace moab; 00063 00064 #define MBI vtkMOABUtils::mbImpl 00065 #define RR if (MB_SUCCESS != result) return result 00066 #define SWAP(a,b) {EntityHandle tmp_ent = a; a = b; b = tmp_ent;} 00067 00068 vtkCellPicker *DrawDual::dualPicker = NULL; 00069 00070 Tag DrawDual::dualCurveTagHandle = 0; 00071 Tag DrawDual::dualSurfaceTagHandle = 0; 00072 00073 DrawDual *DrawDual::gDrawDual = NULL; 00074 Range DrawDual::pickRange; 00075 00076 bool DrawDual::useGraphviz = false; 00077 00078 DrawDual::DrawDual(QLineEdit *pickline1, QLineEdit *pickline2) 00079 : pickLine1(pickline1), pickLine2(pickline2) 00080 { 00081 // make sure we have basic tags we need 00082 EntityHandle dum = 0; 00083 ErrorCode result = MBI->tag_get_handle(DualTool::DUAL_ENTITY_TAG_NAME, 00084 1, MB_TYPE_HANDLE, dualEntityTagHandle, 00085 MB_TAG_DENSE|MB_TAG_CREAT, &dum); 00086 00087 result = MBI->tag_get_handle(DualTool::DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTagHandle); 00088 assert(MB_SUCCESS == result); 00089 00090 result = MBI->tag_get_handle(DualTool::DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTagHandle); 00091 assert(MB_SUCCESS == result); 00092 00093 if (useGraphviz) { 00094 // initialize dot 00095 aginit(); 00096 } 00097 00098 GVEntity *dum2 = NULL; 00099 result = MBI->tag_get_handle("__GVEntity", sizeof(GVEntity*), MB_TAG_OPAQUE, 00100 gvEntityHandle, MB_TAG_DENSE|MB_TAG_CREAT, &dum2); 00101 assert(MB_SUCCESS == result && 0 != gvEntityHandle); 00102 00103 assert(gDrawDual == NULL); 00104 gDrawDual = this; 00105 00106 lastPickedEnt = 0; 00107 secondLastPickedEnt = 0; 00108 00109 add_picker(vtkMOABUtils::myRen); 00110 00111 } 00112 00113 DrawDual::~DrawDual() 00114 { 00115 if (0 != gvEntityHandle) { 00116 Range ents; 00117 MBI->get_entities_by_type_and_tag(0, MBVERTEX, &gvEntityHandle, 00118 NULL, 1, ents, Interface::UNION); 00119 MBI->get_entities_by_type_and_tag(0, MBEDGE, &gvEntityHandle, 00120 NULL, 1, ents, Interface::UNION); 00121 MBI->get_entities_by_type_and_tag(0, MBPOLYGON, &gvEntityHandle, 00122 NULL, 1, ents, Interface::UNION); 00123 std::vector<GVEntity*> gvents(ents.size()); 00124 std::fill(gvents.begin(), gvents.end(), (GVEntity*)NULL); 00125 ErrorCode result = MBI->tag_get_data(gvEntityHandle, ents, &gvents[0]); 00126 if (MB_SUCCESS == result) { 00127 for (unsigned int i = 0; i < gvents.size(); i++) 00128 delete gvents[i]; 00129 } 00130 MBI->tag_delete(gvEntityHandle); 00131 } 00132 00133 if (NULL != dualPicker) { 00134 dualPicker->Delete(); 00135 dualPicker = NULL; 00136 } 00137 00138 if (NULL != vtkMOABUtils::eventCallbackCommand) { 00139 vtkMOABUtils::eventCallbackCommand->Delete(); 00140 vtkMOABUtils::eventCallbackCommand = NULL; 00141 } 00142 00143 if (NULL != gDrawDual) gDrawDual = NULL; 00144 00145 } 00146 00147 void DrawDual::add_picker(vtkRenderer *this_ren) 00148 { 00149 assert(NULL != this_ren); 00150 00151 if (NULL == dualPicker) { 00152 00153 // create a cell picker 00154 dualPicker = vtkCellPicker::New(); 00155 dualPicker->SetTolerance(0.05); 00156 00157 // set up the callback handler for the picker 00158 vtkMOABUtils::eventCallbackCommand = vtkCallbackCommand::New(); 00159 // do we need to do this??? 00160 // eventCallbackCommand->SetClientData(this); 00161 vtkMOABUtils::eventCallbackCommand->SetCallback(process_events); 00162 00163 } 00164 00165 // make sure the renderer can handle observers 00166 vtkRenderWindow *this_ren_window = this_ren->GetRenderWindow(); 00167 assert(NULL != this_ren_window); 00168 vtkRenderWindowInteractor *this_int = this_ren_window->GetInteractor(); 00169 if (NULL == this_int) 00170 this_int = this_ren_window->MakeRenderWindowInteractor(); 00171 00172 assert(NULL != this_int); 00173 00174 vtkInteractorStyle *this_style = vtkInteractorStyle::SafeDownCast( 00175 this_int->GetInteractorStyle()); 00176 assert(NULL != this_style); 00177 this_style->HandleObserversOn(); 00178 00179 // register this command to process the LMB event 00180 this_ren->GetRenderWindow()->GetInteractor()->GetInteractorStyle()-> 00181 AddObserver(vtkCommand::LeftButtonPressEvent, 00182 vtkMOABUtils::eventCallbackCommand); 00183 } 00184 00185 void DrawDual::process_events(vtkObject *caller, 00186 unsigned long event, 00187 void* clientdata, 00188 void* /*vtkNotUsed(calldata)*/) 00189 { 00190 assert(event == vtkCommand::LeftButtonPressEvent); 00191 00192 vtkInteractorStyle* style = reinterpret_cast<vtkInteractorStyle *>(caller); 00193 vtkRenderWindowInteractor *rwi = style->GetInteractor(); 00194 00195 // check to see if the <Ctrl> key is on 00196 if (rwi->GetControlKey()) { 00197 // control key is on - send pick event to picker 00198 if (style->GetState() == VTKIS_NONE) 00199 { 00200 style->FindPokedRenderer(rwi->GetEventPosition()[0], 00201 rwi->GetEventPosition()[1]); 00202 rwi->StartPickCallback(); 00203 vtkRenderer *ren = style->GetCurrentRenderer(); 00204 EntityHandle dual_surf = gDrawDual->get_dual_surf(ren); 00205 if (0 == dual_surf) return; 00206 00207 int ix = rwi->GetEventPosition()[0], 00208 iy = rwi->GetEventPosition()[1]; 00209 00210 double x = (ix - 250) * RAD_PTS/170.0; 00211 double y = iy * RAD_PTS/170.0; 00212 00213 if (my_debug) 00214 std::cout << "Picked point: " << rwi->GetEventPosition()[0] 00215 << ", " << rwi->GetEventPosition()[1] << "(screen), " 00216 << dualPicker->GetMapperPosition()[0] << ", " 00217 << dualPicker->GetMapperPosition()[1] << " (mapper)" 00218 << x << ", " 00219 << y << " (world)" 00220 << std::endl; 00221 00222 Range picked_ents; 00223 ErrorCode result = gDrawDual->process_pick(dual_surf, 00224 x, y, 00225 picked_ents); 00226 if (MB_SUCCESS != result) return; 00227 00228 gDrawDual->print_picked_ents(picked_ents); 00229 00230 if (!picked_ents.empty()) { 00231 // now update the highlighted polydata 00232 } 00233 00234 rwi->EndPickCallback(); 00235 } 00236 } 00237 else { 00238 // otherwise send event to interactor style handler 00239 style->OnLeftButtonDown(); 00240 } 00241 } 00242 00243 void DrawDual::process_pick(vtkRenderer *ren) 00244 { 00245 Range picked_ents; 00246 00247 bool old_picking = false; 00248 if (old_picking) { 00249 00250 assert(0 != dualCurveTagHandle); 00251 00252 // get the actor through the prop, to make sure we get the leaf of an 00253 // assembly 00254 vtkActorCollection *actors = dualPicker->GetActors(); 00255 00256 vtkActor *tmp_actor; 00257 //vtkActor *picked_actor = NULL; 00258 //EntityHandle picked_sheet = 0, picked_chord = 0; 00259 00260 actors->InitTraversal(); 00261 while ((tmp_actor = actors->GetNextItem())) { 00262 EntityHandle this_set = vtkMOABUtils::propSetMap[tmp_actor]; 00263 if (0 == this_set || -1 == dualPicker->GetCellId()) continue; 00264 00265 // get whether it's a dual surface or dual curve 00266 EntityHandle dum_handle = 0; 00267 ErrorCode result = MBI->tag_get_data(dualCurveTagHandle, &this_set, 1, 00268 &dum_handle); 00269 00270 std::cout << "Picked point " << dualPicker->GetCellId() << std::endl; 00271 00272 EntityHandle picked_ent; 00273 if (MB_TAG_NOT_FOUND == result || 0 == dum_handle) 00274 // must be a sheet 00275 picked_ent = gDrawDual->get_picked_cell(this_set, 2, dualPicker->GetCellId()); 00276 else { 00277 picked_ent = gDrawDual->get_picked_cell(this_set, 1, dualPicker->GetCellId()); 00278 } 00279 picked_ents.insert(picked_ent); 00280 } 00281 } 00282 else { 00283 EntityHandle dual_surf = gDrawDual->get_dual_surf(ren); 00284 if (0 == dual_surf) return; 00285 00286 ErrorCode result = gDrawDual->process_pick(dual_surf, 00287 dualPicker->GetPickPosition()[0], 00288 dualPicker->GetPickPosition()[1], 00289 picked_ents); 00290 if (MB_SUCCESS != result) return; 00291 00292 } 00293 00294 if (picked_ents.empty()) return; 00295 00296 gDrawDual->print_picked_ents(picked_ents); 00297 00298 // Range::iterator pit = pickRange.find(picked_ent); 00299 // if (pit == pickRange.end()) pickRange.insert(picked_ent); 00300 // else pickRange.erase(pit); 00301 00302 // now update the highlighted polydata 00303 gDrawDual->update_high_polydatas(); 00304 00305 } 00306 00307 void DrawDual::print_picked_ents(Range &picked_ents, bool from_return) 00308 { 00309 for (Range::iterator rit = picked_ents.begin(); rit != picked_ents.end(); rit++) { 00310 EntityHandle picked_ent = *rit; 00311 00312 // get the vertices 00313 std::ostringstream oss; 00314 const EntityHandle *connect; 00315 int num_connect; 00316 ErrorCode result = MBI->get_connectivity(picked_ent, connect, num_connect); 00317 if (MB_SUCCESS != result) return; 00318 bool first = true; 00319 EntityHandle primals[20]; 00320 std::vector<int> ids; 00321 00322 assert(num_connect < 20); 00323 if (MBI->type_from_handle(picked_ent) == MBPOLYGON) oss << "2-cell: "; 00324 else if (MBI->type_from_handle(picked_ent) == MBEDGE) oss << "1-cell: "; 00325 else oss << "(unknown):"; 00326 result = MBI->tag_get_data(dualEntityTagHandle, connect, num_connect, primals); 00327 ids.resize(num_connect); 00328 result = MBI->tag_get_data(vtkMOABUtils::globalId_tag(), primals, num_connect, &ids[0]); 00329 for (int i = 0; i < num_connect; i++) { 00330 if (!first) oss << "-"; 00331 EntityType this_type = MBI->type_from_handle(primals[i]); 00332 if (this_type == MBHEX) oss << "h"; 00333 else if (this_type == MBQUAD) oss << "f"; 00334 else oss << "u"; 00335 00336 if (ids[i] != 0) oss << ids[i]; 00337 else oss << MBI->id_from_handle(primals[i]); 00338 00339 first = false; 00340 } 00341 00342 if (from_return) { 00343 std::cout << oss.str() << " (" << picked_ent << ")" << std::endl; 00344 pickLine2->setText(QString(oss.str().c_str())); 00345 pickLine1->setText(""); 00346 gDrawDual->secondLastPickedEnt = gDrawDual->lastPickedEnt; 00347 gDrawDual->lastPickedEnt = picked_ent; 00348 00349 gDrawDual->update_high_polydatas(); 00350 } 00351 else if (picked_ent == *picked_ents.begin()) { 00352 //std::cout << oss.str() << " (" << picked_ent << ")" << std::endl; 00353 pickLine2->setText(pickLine1->displayText()); 00354 pickLine1->setText(QString(oss.str().c_str())); 00355 00356 gDrawDual->secondLastPickedEnt = gDrawDual->lastPickedEnt; 00357 gDrawDual->lastPickedEnt = picked_ent; 00358 00359 gDrawDual->update_high_polydatas(); 00360 } 00361 } 00362 } 00363 00364 void DrawDual::update_high_polydatas() 00365 { 00366 // go through each graph window and rebuild picked entities 00367 std::map<EntityHandle, GraphWindows>::iterator mit; 00368 for (mit = surfDrawrings.begin(); mit != surfDrawrings.end(); mit++) { 00369 // reset or initialize 00370 if (NULL == mit->second.pickActor) { 00371 vtkPolyData *pd = vtkPolyData::New(); 00372 vtkPolyDataMapper *mapper = vtkPolyDataMapper::New(); 00373 mapper->SetInput(pd); 00374 mit->second.pickActor = vtkActor::New(); 00375 mit->second.pickActor->SetMapper(mapper); 00376 pd->FastDelete(); 00377 mapper->FastDelete(); 00378 } 00379 else 00380 vtkPolyData::SafeDownCast(mit->second.pickActor->GetMapper()->GetInput())->Reset(); 00381 } 00382 00383 // now go through highlight entities, adding to highlight poly in each window 00384 ErrorCode result; 00385 std::vector<EntityHandle> dual_sets; 00386 std::vector<EntityHandle>::iterator vit; 00387 GVEntity **gvents; 00388 00389 result = MBI->tag_get_data(gvEntityHandle, pickRange, &gvents); 00390 if (MB_SUCCESS != result) return; 00391 unsigned int i, j; 00392 Range::iterator rit; 00393 vtkIdList *id_list = vtkIdList::New(); 00394 id_list->Allocate(1); 00395 00396 for (i = 0, rit = pickRange.begin(); rit != pickRange.end(); i++, rit++) { 00397 // can be up to 3 instances of this entity 00398 for (j = 0; j < 3; j++) { 00399 if (gvents[i]->myActors[j] == NULL) continue; 00400 00401 // get the vtk cell and duplicate it 00402 vtkPolyData *this_pd = 00403 vtkPolyData::SafeDownCast(gvents[i]->myActors[j]->GetMapper()->GetInput()); 00404 assert(NULL != this_pd); 00405 assert(gvents[i]->dualSurfs[j] != 0); 00406 GraphWindows &gw = surfDrawrings[gvents[i]->dualSurfs[j]]; 00407 vtkPolyData *pd = vtkPolyData::SafeDownCast(gw.pickActor->GetMapper()->GetInput()); 00408 pd->Allocate(); 00409 id_list->Reset(); 00410 id_list->InsertNextId(gvents[i]->vtkEntityIds[j]); 00411 pd->CopyCells(this_pd, id_list); 00412 00413 00414 vtkInteractorStyle *this_style = vtkInteractorStyle::SafeDownCast( 00415 gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->GetInteractor()->GetInteractorStyle()); 00416 this_style->HighlightProp(gw.pickActor); 00417 } 00418 } 00419 id_list->FastDelete(); 00420 } 00421 00422 EntityHandle DrawDual::get_picked_cell(EntityHandle cell_set, 00423 const int dim, 00424 const int picked_cell) 00425 { 00426 // get the cells in the set 00427 Range cells; 00428 ErrorCode result; 00429 if (1 == dim) 00430 result = vtkMOABUtils::dualTool->get_dual_entities(cell_set, NULL, &cells, NULL, NULL, NULL); 00431 else if (2 == dim) 00432 result = vtkMOABUtils::dualTool->get_dual_entities(cell_set, &cells, NULL, NULL, NULL, NULL); 00433 else 00434 assert(false); 00435 00436 if (MB_SUCCESS != result) return 0; 00437 00438 Range::iterator rit = cells.begin(); 00439 00440 for (int i = 0; i < picked_cell; i++, rit++); 00441 00442 return *rit; 00443 } 00444 00445 bool DrawDual::print_dual_surfs(Range &dual_surfs, 00446 const bool /*use_offsets*/) 00447 { 00448 ErrorCode success = MB_SUCCESS; 00449 int offset = 0; 00450 00451 // create vtkWindowToImageFilter 00452 vtkWindowToImageFilter *wtif = vtkWindowToImageFilter::New(); 00453 00454 // create vtkPNGWriter 00455 vtkPNGWriter *pngw = vtkPNGWriter::New(); 00456 00457 // set vtkPNGWriter input to output port of vtkWindowToImageFilter 00458 pngw->SetInput(wtif->GetOutput()); 00459 00460 char filename[15]; 00461 std::vector<int> surf_ids(dual_surfs.size()); 00462 success = vtkMOABUtils::MBI->tag_get_data(vtkMOABUtils::globalId_tag(), 00463 dual_surfs, &surf_ids[0]); 00464 if (MB_SUCCESS != success) return success; 00465 Range::iterator rit; 00466 int i; 00467 00468 for (rit = dual_surfs.begin(), i = 0; 00469 rit != dual_surfs.end(); rit++, i++) { 00470 00471 GraphWindows &this_gw = surfDrawrings[*rit]; 00472 00473 EntityHandle dum_handle = 0; 00474 ErrorCode tmp_success = MBI->tag_get_data(dualSurfaceTagHandle, &(*rit), 1, 00475 &dum_handle); 00476 if (MB_TAG_NOT_FOUND == tmp_success || dum_handle == 0) continue; 00477 00478 tmp_success = draw_dual_surf(*rit, offset); 00479 if (MB_SUCCESS != tmp_success) success = tmp_success; 00480 // if (use_offsets) offset++; 00481 00482 // set vtkWindowToImageFilter input to render window 00483 wtif->SetInput(this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()); 00484 wtif->Update(); 00485 00486 // set file name 00487 sprintf(filename, "s%d.png", surf_ids[i]); 00488 pngw->SetFileName(filename); 00489 00490 // call Write function on png writer 00491 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->Render(); 00492 pngw->Write(); 00493 } 00494 00495 wtif->Delete(); 00496 pngw->Delete(); 00497 00498 return (MB_SUCCESS == success ? true : false); 00499 } 00500 00501 bool DrawDual::draw_dual_surfs(Range &dual_surfs, 00502 const bool /*use_offsets*/) 00503 { 00504 ErrorCode success = MB_SUCCESS; 00505 int offset = 0; 00506 for (Range::reverse_iterator rit = dual_surfs.rbegin(); rit != dual_surfs.rend(); rit++) { 00507 EntityHandle dum_handle = 0; 00508 ErrorCode tmp_success = MBI->tag_get_data(dualSurfaceTagHandle, &(*rit), 1, 00509 &dum_handle); 00510 if (MB_TAG_NOT_FOUND == tmp_success || dum_handle == 0) continue; 00511 00512 tmp_success = draw_dual_surf(*rit, offset); 00513 if (MB_SUCCESS != tmp_success) success = tmp_success; 00514 // if (use_offsets) offset++; 00515 } 00516 00517 return (MB_SUCCESS == success ? true : false); 00518 } 00519 00520 bool DrawDual::draw_dual_surfs(std::vector<EntityHandle> &dual_surfs, 00521 const bool /*use_offsets*/) 00522 { 00523 ErrorCode success = MB_SUCCESS; 00524 int offset = 0; 00525 for (std::vector<EntityHandle>::reverse_iterator vit = dual_surfs.rbegin(); 00526 vit != dual_surfs.rend(); vit++) { 00527 EntityHandle dum_handle = 0; 00528 ErrorCode tmp_success = MBI->tag_get_data(dualSurfaceTagHandle, &(*vit), 1, 00529 &dum_handle); 00530 if (MB_TAG_NOT_FOUND == tmp_success || dum_handle == 0) continue; 00531 00532 tmp_success = draw_dual_surf(*vit, offset); 00533 if (MB_SUCCESS != tmp_success) success = tmp_success; 00534 //if (use_offsets) offset++; 00535 } 00536 00537 return (MB_SUCCESS == success ? true : false); 00538 } 00539 00540 ErrorCode DrawDual::draw_dual_surf(EntityHandle dual_surf, 00541 int offset_num) 00542 { 00543 // draw this dual surface 00544 00545 // 3. make a new pd for this drawring 00546 GraphWindows &this_gw = surfDrawrings[dual_surf]; 00547 vtkPolyData *pd; 00548 get_clean_pd(dual_surf, this_gw.sheetDiagram, pd); 00549 00550 this_gw.sheetDiagram->show(); 00551 00552 // get the cells and vertices on this dual surface 00553 Range dcells, dedges, dverts, face_verts, loop_edges; 00554 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, &dcells, &dedges, 00555 &dverts, &face_verts, &loop_edges); 00556 if (MB_SUCCESS != result) return result; 00557 00558 if (dcells.empty() || dedges.empty() || dverts.empty()) return MB_FAILURE; 00559 00560 // 1. gather/construct data for graphviz 00561 ErrorCode success = construct_graphviz_data(dual_surf, dcells, dedges, dverts, 00562 face_verts, loop_edges); 00563 if (MB_SUCCESS != success) return success; 00564 00565 00566 // 2. tell graphviz to compute graph 00567 // gvBindContext((GVC_t*)gvizGvc, this_gw.gvizGraph); 00568 // neato_layout(this_gw.gvizGraph); 00569 if (my_debug) { 00570 std::cout << "Before layout:" << std::endl; 00571 if (useGraphviz) agwrite(this_gw.gvizGraph, stdout); 00572 } 00573 // neato_init_graph(this_gw.gvizGraph); 00574 if (useGraphviz) gvLayout(gvContext(), this_gw.gvizGraph, (char*)"neato"); 00575 else smooth_dual_surf(dual_surf, dcells, dedges, dverts, 00576 face_verts, loop_edges); 00577 00578 if (my_debug) { 00579 std::cout << "After layout, before vtk:" << std::endl; 00580 if (useGraphviz) agwrite(this_gw.gvizGraph, stdout); 00581 } 00582 00583 success = fixup_degen_bchords(dual_surf); 00584 if (MB_SUCCESS != success) return success; 00585 00586 vtkRenderer *this_ren = 00587 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->GetRenderers()->GetFirstRenderer(); 00588 00589 // 4. create vtk points, cells for this 2d dual surface drawing 00590 success = make_vtk_data(dual_surf, pd, this_ren); 00591 if (MB_SUCCESS != success) return success; 00592 if (my_debug) { 00593 std::cout << "After layout, after vtk:" << std::endl; 00594 if (useGraphviz) agwrite(this_gw.gvizGraph, stdout); 00595 } 00596 00597 // 5. generate "other sheet" labels 00598 vtkPolyData *new_pd; 00599 result = label_other_sheets(dual_surf, pd, new_pd); 00600 if (MB_SUCCESS != result) return result; 00601 00602 pd->Update(); 00603 00604 //this_gw.qvtkWidget->GetRenderWindow()->DebugOn(); 00605 00606 // 5. render the window 00607 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->Render(); 00608 00609 // now that we've rendered, can get the size of various things 00610 int *tmp_siz = this_ren->GetSize(); 00611 xSize = tmp_siz[0]; 00612 ySize = tmp_siz[1]; 00613 tmp_siz = this_ren->GetOrigin(); 00614 xOrigin = tmp_siz[0]; 00615 yOrigin = tmp_siz[1]; 00616 00617 // 6. draw labels for dual surface, points, dual curves 00618 success = draw_labels(dual_surf, pd, new_pd); 00619 if (MB_SUCCESS != success) return success; 00620 pd->FastDelete(); 00621 new_pd->FastDelete(); 00622 00623 // 7. add a picker 00624 add_picker(this_ren); 00625 00626 if (my_debug && useGraphviz) agwrite(this_gw.gvizGraph, stdout); 00627 00628 int old_pos[2] = {0, 0}; 00629 00630 int *win_siz = this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->GetSize(); 00631 int new_pos[2] = {(int) ((double)old_pos[0] + win_siz[0]*(.1*offset_num)), old_pos[1]}; 00632 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->SetPosition(new_pos); 00633 00634 return success; 00635 } 00636 00637 ErrorCode DrawDual::fixup_degen_bchords(EntityHandle dual_surf) 00638 { 00639 // make sure the mid-pt of degenerate blind chord segments is on the correct 00640 // side of the other chord it splits 00641 Range dcells, dedges, dverts, face_verts, loop_edges; 00642 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, &dcells, &dedges, 00643 &dverts, &face_verts, &loop_edges); RR; 00644 00645 Range tmp_edges, degen_2cells; 00646 00647 double avg_pos0[3], avg_pos1[3], dum_pos0[3], dum_pos1[3], dum_pos2[3]; 00648 00649 for (Range::iterator rit = dcells.begin(); rit != dcells.end(); rit++) { 00650 // first, find if it's degenerate 00651 tmp_edges.clear(); 00652 result = MBI->get_adjacencies(&(*rit), 1, 1, false, tmp_edges); RR; 00653 if (tmp_edges.size() != 2) continue; 00654 // also skip if it's on the boundary 00655 else if (loop_edges.find(*tmp_edges.begin()) != loop_edges.end() || 00656 loop_edges.find(*tmp_edges.rbegin()) != loop_edges.end()) 00657 continue; 00658 00659 degen_2cells.insert(*rit); 00660 } 00661 00662 MeshTopoUtil mtu(vtkMOABUtils::mbImpl); 00663 00664 while (!degen_2cells.empty()) { 00665 // grab the first one off the list 00666 EntityHandle tcell1 = *degen_2cells.begin(); 00667 00668 // look for any adjacent degenerate 2cells also on the list 00669 const EntityHandle *connect; 00670 int num_connect; 00671 Range adj_2cells; 00672 result = MBI->get_connectivity(tcell1, connect, num_connect); RR; 00673 result = MBI->get_adjacencies(connect, num_connect, 2, false, adj_2cells); RR; 00674 adj_2cells = intersect( degen_2cells, adj_2cells); 00675 if (!adj_2cells.empty()) degen_2cells = subtract( degen_2cells, adj_2cells); 00676 00677 // ok, have all the adjacent degen 2cells; get the 1cells 00678 Range adj_1cells; 00679 result = vtkMOABUtils::mbImpl->get_adjacencies(adj_2cells, 1, false, adj_1cells, 00680 Interface::UNION); 00681 if (MB_SUCCESS != result) return result; 00682 00683 // decide whether this sheet is blind or not 00684 bool sheet_is_blind = vtkMOABUtils::dualTool->is_blind(dual_surf); 00685 00686 // branch depending on what kind of arrangement we have 00687 if (adj_2cells.size() == 2 && !sheet_is_blind) { 00688 assert(3 == adj_1cells.size()); 00689 // blind chord intersecting another chord 00690 // get the middle edge and chord, and the next 2 verts along the chord 00691 Range dum; 00692 result = MBI->get_adjacencies(adj_2cells, 1, false, dum); RR; 00693 assert(1 == dum.size()); 00694 EntityHandle middle_edge = *dum.begin(); 00695 EntityHandle chord = vtkMOABUtils::dualTool->get_dual_hyperplane(middle_edge); 00696 assert(0 != chord); 00697 EntityHandle verts[2]; 00698 result = vtkMOABUtils::dualTool->get_opposite_verts(middle_edge, chord, verts); RR; 00699 00700 // get the gv points for the four vertices 00701 void *next_points[2], *points[2]; 00702 get_graph_points(verts, 2, dual_surf, next_points); 00703 assert(next_points[0] != NULL && next_points[1] != NULL); 00704 result = MBI->get_connectivity(middle_edge, connect, num_connect); RR; 00705 get_graph_points(connect, 2, dual_surf, points); 00706 assert(points[0] != NULL && points[1] != NULL); 00707 00708 // now space points along the line segment joining the next_points 00709 get_graphpoint_pos(next_points[0], dum_pos0); 00710 get_graphpoint_pos(next_points[1], dum_pos1); 00711 00712 avg_pos0[0] = .5*(dum_pos0[0] + dum_pos1[0]); 00713 avg_pos0[1] = .5*(dum_pos0[1] + dum_pos1[1]); 00714 dum_pos0[0] = 0.5 * (avg_pos0[0] + dum_pos0[0]); 00715 dum_pos0[1] = 0.5 * (avg_pos0[1] + dum_pos0[1]); 00716 dum_pos1[0] = 0.5 * (avg_pos0[0] + dum_pos1[0]); 00717 dum_pos1[1] = 0.5 * (avg_pos0[1] + dum_pos1[1]); 00718 set_graphpoint_pos(points[0], dum_pos0); 00719 set_graphpoint_pos(points[1], dum_pos1); 00720 00721 // also fix the middle point on this edge 00722 points[0] = NULL; 00723 get_graph_points(&middle_edge, 1, dual_surf, points); 00724 set_graphpoint_pos(points[0], avg_pos0); 00725 00726 // now fix the other 2 dedges 00727 adj_1cells.erase(middle_edge); 00728 for (Range::iterator rit = adj_1cells.begin(); rit != adj_1cells.end(); rit++) { 00729 // get the other 2cell 00730 Range dum = dcells; 00731 result = MBI->get_adjacencies(&(*rit), 1, 2, false, dum); 00732 dum = subtract( dum, adj_2cells); 00733 assert(1 == dum.size()); 00734 // get the vertices and points of them, and average their positions 00735 const EntityHandle *connect2; 00736 result = MBI->get_connectivity(*dum.begin(), connect2, num_connect); RR; 00737 std::vector<void*> tc_points(num_connect); 00738 get_graph_points(connect2, num_connect, dual_surf, &tc_points[0]); 00739 avg_pos1[0] = avg_pos1[1] = avg_pos1[2] = 0.0; 00740 for (int i = 0; i < num_connect; i++) { 00741 if (connect2[i] != connect[0] && connect2[i] != connect[1]) { 00742 get_graphpoint_pos(tc_points[i], dum_pos0); 00743 avg_pos1[0] += dum_pos0[0]; 00744 avg_pos1[1] += dum_pos0[1]; 00745 } 00746 } 00747 avg_pos1[0] = (.4*avg_pos1[0]/(num_connect-2) + .6*avg_pos0[0]); 00748 avg_pos1[1] = (.4*avg_pos1[1]/(num_connect-2) + .6*avg_pos0[1]); 00749 get_graph_points(&(*rit), 1, dual_surf, &tc_points[0]); 00750 set_graphpoint_pos(tc_points[0], avg_pos1); 00751 } 00752 } 00753 else if (adj_2cells.size() == 1) { 00754 assert(adj_1cells.size() == 2); 00755 00756 // get vertices making up degen 2cell and their avg position 00757 const EntityHandle *connect; 00758 result = MBI->get_connectivity(*adj_2cells.begin(), connect, num_connect); RR; 00759 std::vector<void*> tc_points(num_connect); 00760 get_graph_points(connect, num_connect, dual_surf, &tc_points[0]); 00761 get_graphpoint_pos(tc_points[0], dum_pos0); 00762 get_graphpoint_pos(tc_points[1], dum_pos1); 00763 00764 avg_pos0[0] = .5*(dum_pos0[0] + dum_pos1[0]); 00765 avg_pos0[1] = .5*(dum_pos0[1] + dum_pos1[1]); 00766 00767 // for each 1cell, get the vertices on the adjacent non-degen 2cell 00768 // and points of them, and average their positions 00769 for (Range::iterator rit = adj_1cells.begin(); rit != adj_1cells.end(); rit++) { 00770 // get the other 2cell 00771 Range dum = dcells; 00772 result = MBI->get_adjacencies(&(*rit), 1, 2, false, dum); 00773 dum = subtract( dum, adj_2cells); 00774 assert(1 == dum.size()); 00775 // get the vertices and points of them, and average their positions 00776 const EntityHandle *connect; 00777 result = MBI->get_connectivity(*dum.begin(), connect, num_connect); RR; 00778 std::vector<void*> tc_points(num_connect); 00779 get_graph_points(connect, num_connect, dual_surf, &tc_points[0]); 00780 avg_pos1[0] = avg_pos1[1] = avg_pos1[2] = 0.0; 00781 for (int i = 0; i < num_connect; i++) { 00782 get_graphpoint_pos(tc_points[i], dum_pos0); 00783 avg_pos1[0] += dum_pos0[0]; 00784 avg_pos1[1] += dum_pos0[1]; 00785 } 00786 avg_pos1[0] = (.2*avg_pos1[0]/num_connect + .8*avg_pos0[0]); 00787 avg_pos1[1] = (.2*avg_pos1[1]/num_connect + .8*avg_pos0[1]); 00788 get_graph_points(&(*rit), 1, dual_surf, &tc_points[0]); 00789 set_graphpoint_pos(tc_points[0], avg_pos1); 00790 } 00791 } 00792 else if ((adj_2cells.size() == 4 && adj_1cells.size() == 4) || 00793 sheet_is_blind) { 00794 // pillow sheet, right after atomic pillow; just place 1cell mid-pts so 00795 // we can see them 00796 // get # chords, since drawing method depends on that 00797 Range chords; 00798 result = MBI->get_child_meshsets(dual_surf, chords); 00799 if (MB_SUCCESS != result) return result; 00800 00801 if (2 == chords.size()) { 00802 00803 const EntityHandle *connect; 00804 int num_connect; 00805 result = MBI->get_connectivity(*adj_1cells.begin(), connect, num_connect); RR; 00806 void *vert_pts[2], *edge_pts[4]; 00807 get_graph_points(connect, 2, dual_surf, vert_pts); 00808 std::vector<EntityHandle> edges; 00809 std::copy(adj_1cells.begin(), adj_1cells.end(), std::back_inserter(edges)); 00810 00811 // check that 1cells are in proper order, i.e. they share 2cell with adj 1cell 00812 Range dum; 00813 result = MBI->get_adjacencies(&edges[0], 2, 2, false, dum); RR; 00814 dum = intersect( dum, adj_2cells); 00815 if (dum.empty()) { 00816 // not adjacent - switch list order 00817 EntityHandle dum_h = edges[1]; 00818 edges[1] = edges[2]; 00819 edges[2] = dum_h; 00820 } 00821 00822 get_graph_points(&edges[0], 4, dual_surf, edge_pts); 00823 dum_pos0[0] = CENT_X; dum_pos0[1] = CENT_Y+RAD_PTS; 00824 dum_pos1[0] = CENT_X; dum_pos1[1] = CENT_Y-RAD_PTS; 00825 set_graphpoint_pos(vert_pts[0], dum_pos0); 00826 set_graphpoint_pos(vert_pts[1], dum_pos1); 00827 for (int i = 0; i < 4; i++) { 00828 dum_pos0[1] = CENT_Y; dum_pos0[0] = CENT_X + (2*i-3)*2*RAD_PTS/5.0; 00829 set_graphpoint_pos(edge_pts[i], dum_pos0); 00830 } 00831 } 00832 else if (3 == chords.size()) { 00833 // get the middle chord 00834 EntityHandle middle_chord = 0; 00835 for (Range::iterator rit = chords.begin(); rit != chords.end(); rit++) { 00836 int num_ents; 00837 result = MBI->get_number_entities_by_type(*rit, MBEDGE, num_ents); 00838 if (MB_SUCCESS != result) return result; 00839 if (2 < num_ents) { 00840 middle_chord = *rit; 00841 break; 00842 } 00843 } 00844 if (0 == middle_chord) return MB_FAILURE; 00845 00846 chords.erase(middle_chord); 00847 00848 // get the edges on each of the non-middle chords 00849 std::vector<EntityHandle> chord_edges[2]; 00850 result = MBI->get_entities_by_handle(*chords.begin(), chord_edges[0]); RR; 00851 result = MBI->get_entities_by_handle(*chords.rbegin(), chord_edges[1]); RR; 00852 00853 // align them such that chord_edges[0][0] and chord_edges[1][0] do not share a 2cell 00854 // on this sheet; arbitrarily choose chord_edges[0][0] to be left-most 00855 EntityHandle shared_2cell = mtu.common_entity(chord_edges[0][0], chord_edges[1][0], 2); 00856 if (0 != shared_2cell && vtkMOABUtils::dualTool->get_dual_hyperplane(shared_2cell) == dual_surf) { 00857 shared_2cell = chord_edges[0][0]; 00858 chord_edges[0][0] = chord_edges[0][1]; 00859 chord_edges[0][1] = shared_2cell; 00860 } 00861 if (0 != mtu.common_entity(chord_edges[0][0], chord_edges[1][0], 2)) 00862 return MB_FAILURE; 00863 00864 double num_x = 2; 00865 double xdelta = (RAD_PTS/num_x)/3.0, xcent = CENT_X; 00866 double xpos = CENT_X - xdelta; 00867 00868 for (int i = 0; i < num_x; i++) { 00869 00870 // get the edge on the middle chord between chord_edges[i][0] and chord_edges[i][1]; that will 00871 // be the intersection of edges on middle chord and edges adjacent to vertices 00872 // bounding chord_edges[i][0] 00873 Range middle_edges; 00874 result = MBI->get_entities_by_handle(middle_chord, middle_edges); RR; 00875 const EntityHandle *connect; 00876 int num_connect; 00877 result = MBI->get_connectivity(*chord_edges[i].begin(), connect, num_connect); RR; 00878 result = MBI->get_adjacencies(connect, num_connect, 1, false, middle_edges); RR; 00879 if (1 != middle_edges.size()) { 00880 // take off any edges that aren't on this sheet? no, just punt... 00881 std::cout << "Warning: not a clear middle edge..." << std::endl; 00882 } 00883 00884 // get the points for the two vertices and the 3 edges 00885 // non-middle chord; get the edges too 00886 void *vert_pts[2], *edge_pts[3]; 00887 get_graph_points(connect, 2, dual_surf, vert_pts); 00888 get_graph_points(&chord_edges[i][0], 2, dual_surf, edge_pts); 00889 get_graph_points(&(*middle_edges.begin()), 1, dual_surf, &edge_pts[2]); 00890 00891 dum_pos0[0] = xpos; xpos += xdelta; 00892 dum_pos2[0] = (xpos < xcent ? xpos-xdelta/2 : xpos+xdelta/2); 00893 avg_pos0[0] = xpos; 00894 avg_pos1[0] = xpos; xpos += xdelta; 00895 dum_pos1[0] = xpos; xpos += xdelta; 00896 00897 avg_pos0[1] = CENT_Y - .5*RAD_PTS; 00898 avg_pos1[1] = CENT_Y + .5*RAD_PTS; 00899 dum_pos0[1] = dum_pos1[1] = dum_pos2[1] = CENT_Y; 00900 00901 set_graphpoint_pos(vert_pts[0], avg_pos0); 00902 set_graphpoint_pos(vert_pts[1], avg_pos1); 00903 set_graphpoint_pos(edge_pts[0], dum_pos0); 00904 set_graphpoint_pos(edge_pts[1], dum_pos1); 00905 set_graphpoint_pos(edge_pts[2], dum_pos2); 00906 00907 xpos += xdelta; 00908 } 00909 } 00910 } 00911 } 00912 00913 return MB_SUCCESS; 00914 } 00915 00916 ErrorCode DrawDual::get_graphpoint_pos(void *point, double *pos) 00917 { 00918 if (useGraphviz) { 00919 Agnode_t *this_node = (Agnode_t*) point; 00920 pos[0] = ND_coord_i(this_node).x; 00921 pos[1] = ND_coord_i(this_node).y; 00922 pos[2] = 0.0; 00923 } 00924 else { 00925 EntityHandle this_node = (EntityHandle) point; 00926 ErrorCode result = MBI->get_coords(&this_node, 1, pos); 00927 if (MB_SUCCESS != result) return result; 00928 } 00929 00930 return MB_SUCCESS; 00931 } 00932 00933 ErrorCode DrawDual::set_graphpoint_pos(void *point, double *pos) 00934 { 00935 if (useGraphviz) { 00936 Agnode_t *this_node = (Agnode_t*) point; 00937 ND_coord_i(this_node).x = pos[0]; 00938 ND_coord_i(this_node).y = pos[1]; 00939 } 00940 else { 00941 pos[2] = 0.0; 00942 EntityHandle this_node = (EntityHandle) point; 00943 ErrorCode result = MBI->set_coords(&this_node, 1, pos); 00944 if (MB_SUCCESS != result) return result; 00945 } 00946 00947 return MB_SUCCESS; 00948 } 00949 00950 ErrorCode DrawDual::make_vtk_data(EntityHandle dual_surf, 00951 vtkPolyData *pd, 00952 vtkRenderer *this_ren) 00953 { 00954 // get the cells and vertices on this dual surface 00955 Range dcells, dverts, dverts_loop, dedges, dedges_loop; 00956 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, &dcells, &dedges, 00957 &dverts, &dverts_loop, &dedges_loop); 00958 if (MB_SUCCESS != result) return result; 00959 00960 // get the GVEntity for each entity 00961 std::map<EntityHandle, GVEntity *> vert_gv_map; 00962 00963 // allocate cells and points; start by getting the point and cell arrays 00964 result = allocate_points(dual_surf, pd, dverts, dedges, vert_gv_map); 00965 if (MB_SUCCESS != result) return result; 00966 00967 // get an int field to put the color id into 00968 vtkFloatArray *color_ids = vtkFloatArray::New(); 00969 color_ids->SetName("ColorId"); 00970 color_ids->Initialize(); 00971 00972 // make the 2d cells 00973 int global_id; 00974 pd->Allocate(); 00975 result = vtkMOABUtils::MBI->tag_get_data(vtkMOABUtils::globalId_tag(), &dual_surf, 00976 1, &global_id); 00977 if (MB_SUCCESS != result) return result; 00978 00979 result = make_vtk_cells(dcells, 2, (float) global_id, 00980 dual_surf, vert_gv_map, pd, 00981 color_ids); 00982 if (MB_SUCCESS != result) return result; 00983 00984 // make the color ids the active scalar 00985 pd->GetCellData()->AddArray(color_ids); 00986 pd->GetCellData()->SetScalars(color_ids); 00987 pd->GetCellData()->SetActiveAttribute("ColorId", 0); 00988 color_ids->FastDelete(); 00989 00990 // make the 1d cells chord by chord 00991 std::vector<EntityHandle> chords; 00992 result = MBI->get_child_meshsets(dual_surf, chords); 00993 if (MB_SUCCESS != result) return result; 00994 00995 for (std::vector<EntityHandle>::iterator vit = chords.begin(); 00996 vit != chords.end(); vit++) { 00997 // set color of chord to other sheet's color 00998 EntityHandle color_set = other_sheet(*vit, dual_surf); 00999 result = vtkMOABUtils::MBI->tag_get_data(vtkMOABUtils::globalId_tag(), &color_set, 01000 1, &global_id); 01001 if (MB_SUCCESS != result) return result; 01002 01003 // get edges in this chord 01004 Range dedges; 01005 result = vtkMOABUtils::dualTool->get_dual_entities(*vit, NULL, &dedges, NULL, NULL, NULL); 01006 if (MB_SUCCESS != result) return result; 01007 01008 // construct a new polydata, and borrow the points from the sheet's pd 01009 vtkPolyData *chord_pd = vtkPolyData::New(); 01010 chord_pd->Allocate(); 01011 chord_pd->SetPoints(pd->GetPoints()); 01012 vtkPolyDataMapper *chord_mapper = vtkPolyDataMapper::New(); 01013 chord_mapper->SetInput(chord_pd); 01014 vtkActor *chord_actor = vtkActor::New(); 01015 vtkMOABUtils::propSetMap[chord_actor] = *vit; 01016 chord_actor->SetMapper(chord_mapper); 01017 this_ren->AddActor(chord_actor); 01018 double red, green, blue; 01019 vtkMOABUtils::get_colors(color_set, 0, global_id, red, green, blue); 01020 chord_actor->GetProperty()->SetColor(red, green, blue); 01021 chord_actor->GetProperty()->SetLineWidth(2.0); 01022 01023 // now make the 1-cells 01024 result = make_vtk_cells(dedges, 1, (float) global_id, 01025 dual_surf, vert_gv_map, chord_pd, NULL); 01026 chord_pd->FastDelete(); 01027 chord_actor->FastDelete(); 01028 chord_mapper->FastDelete(); 01029 if (MB_SUCCESS != result) return result; 01030 } 01031 01032 return MB_SUCCESS; 01033 } 01034 01035 ErrorCode DrawDual::make_vtk_cells(const Range &cell_range, const int dim, 01036 const float color_index, 01037 const EntityHandle dual_surf, 01038 std::map<EntityHandle, GVEntity *> &vert_gv_map, 01039 vtkPolyData *pd, 01040 vtkFloatArray *color_ids) 01041 { 01042 std::vector<EntityHandle> cell_verts; 01043 std::vector<GVEntity *> gv_cells(cell_range.size()); 01044 int cell_num; 01045 Range::iterator rit; 01046 vtkIdType cell_points[20]; 01047 static int vtk_cell_type[] = {VTK_VERTEX, VTK_LINE, VTK_POLYGON}; 01048 01049 ErrorCode result = MBI->tag_get_data(gvEntityHandle, cell_range, &gv_cells[0]); 01050 if (MB_SUCCESS != result) return result; 01051 01052 Range cell_edges, shared_edges, tmp_verts; 01053 for (rit = cell_range.begin(), cell_num = 0; 01054 rit != cell_range.end(); rit++, cell_num++) { 01055 // get the vertices in this cell; must be done through vector for polygons 01056 cell_verts.clear(); 01057 result = MBI->get_adjacencies(&(*rit), 1, 0, false, cell_verts); RR; 01058 assert(cell_verts.size() <= 20); 01059 cell_edges.clear(); 01060 result = MBI->get_adjacencies(&(*rit), 1, 1, false, cell_edges); RR; 01061 01062 // for each, check for vtk point & build one if it's missing, 01063 // then create the vtk entity 01064 01065 // in this loop, num_pts keeps the place in the edge/2cell's vtk point 01066 // vector (cell_points), while i keeps the place in the (linear) vertex list 01067 // for the edge/2cell 01068 int num_pts = 0; 01069 for (unsigned int i = 0; i < cell_verts.size(); i++) { 01070 GVEntity *this_gv = vert_gv_map[cell_verts[i]]; 01071 assert(this_gv != NULL); 01072 int index = this_gv->get_index(dual_surf); 01073 assert(index >= 0 && NULL != this_gv->gvizPoints[index]); 01074 cell_points[num_pts++] = this_gv->vtkEntityIds[index]; 01075 01076 // check for higher-order edge 01077 shared_edges.clear(); 01078 tmp_verts.clear(); 01079 tmp_verts.insert(cell_verts[i]); 01080 tmp_verts.insert(cell_verts[(i+1)%cell_verts.size()]); 01081 result = MBI->get_adjacencies(tmp_verts, 1, false, shared_edges); RR; 01082 if (shared_edges.size() > 1) { 01083 // filter for ones in this cell 01084 shared_edges = intersect( shared_edges, cell_edges); 01085 assert(!shared_edges.empty()); 01086 // get the mid-pt of this edge and include in list; if we're inside a 2-edge 01087 // cell and we're on the 2nd vertex, take the 2nd edge 01088 EntityHandle this_edge = *shared_edges.begin(); 01089 if (cell_verts.size() == 2 && i != 0) this_edge = *shared_edges.rbegin(); 01090 this_gv = vert_gv_map[this_edge]; 01091 assert(this_gv != NULL); 01092 index = this_gv->get_index(dual_surf); 01093 assert(index >= 0 && this_gv->gvizPoints[index] != NULL); 01094 cell_points[num_pts++] = this_gv->vtkEntityIds[index+2]; 01095 } 01096 } 01097 01098 if (dim == 2) 01099 cell_points[num_pts++] = cell_points[0]; 01100 01101 // create the vtk cell 01102 if (NULL == gv_cells[cell_num]) { 01103 gv_cells[cell_num] = new GVEntity(); 01104 gv_cells[cell_num]->moabEntity = *rit; 01105 } 01106 01107 int index = gv_cells[cell_num]->get_index(dual_surf); 01108 assert(index > -10); 01109 if (index < 0) { 01110 index = -index - 1; 01111 gv_cells[cell_num]->dualSurfs[index] = dual_surf; 01112 } 01113 01114 int this_type = vtk_cell_type[dim]; 01115 if (dim == 1 && num_pts > 2) this_type = VTK_POLY_LINE; 01116 int this_cell = pd->InsertNextCell(this_type, num_pts, 01117 cell_points); 01118 gv_cells[cell_num]->vtkEntityIds[index] = this_cell; 01119 if (NULL != color_ids) 01120 color_ids->InsertValue(this_cell, color_index); 01121 } 01122 01123 result = MBI->tag_set_data(gvEntityHandle, cell_range, &gv_cells[0]); 01124 if (MB_SUCCESS != result) return result; 01125 01126 return MB_SUCCESS; 01127 } 01128 01129 ErrorCode DrawDual::allocate_points(EntityHandle dual_surf, 01130 vtkPolyData *pd, 01131 Range &verts, 01132 Range &edges, 01133 std::map<EntityHandle, GVEntity*> &vert_gv_map) 01134 { 01135 vtkPoints *points = pd->GetPoints(); 01136 if (NULL == points) { 01137 points = vtkPoints::New(); 01138 pd->SetPoints(points); 01139 points->Delete(); 01140 } 01141 pd->Allocate(); 01142 01143 std::vector<GVEntity*> gv_verts; 01144 gv_verts.resize(verts.size()); 01145 01146 // get the gventities 01147 ErrorCode result = MBI->tag_get_data(gvEntityHandle, verts, &gv_verts[0]); 01148 if (MB_SUCCESS != result) { 01149 std::cout << "Couldn't get gv vertex data." << std::endl; 01150 return result; 01151 } 01152 01153 unsigned int i; 01154 Range::iterator rit; 01155 char dum_str[80]; 01156 Agsym_t *asym_pos = (useGraphviz ? get_asym(dual_surf, 0, "pos") : NULL); 01157 double dum_pos[3]; 01158 01159 // get the transform based on old/new positions for loop vert(s) 01160 double x_xform = 1.0, y_xform = 1.0; 01161 if (useGraphviz) 01162 result = get_xform(dual_surf, asym_pos, x_xform, y_xform); 01163 01164 for (rit = verts.begin(), i = 0; rit != verts.end(); rit++, i++) { 01165 assert(NULL != gv_verts[i]); 01166 int index = gv_verts[i]->get_index(dual_surf); 01167 assert(index >= 0 && NULL != gv_verts[i]->gvizPoints[index]); 01168 if (gv_verts[i]->vtkEntityIds[index] == -1) { 01169 get_graphpoint_pos(gv_verts[i]->gvizPoints[index], dum_pos); 01170 if (useGraphviz) { 01171 sprintf(dum_str, " %d,%d", (int)dum_pos[0], (int)dum_pos[1]); 01172 agxset(gv_verts[i]->gvizPoints[index], asym_pos->index, dum_str); 01173 } 01174 01175 gv_verts[i]->vtkEntityIds[index] = 01176 points->InsertNextPoint(dum_pos[0], dum_pos[1], dum_pos[2]); 01177 } 01178 vert_gv_map[*rit] = gv_verts[i]; 01179 } 01180 01181 if (edges.empty()) return MB_SUCCESS; 01182 01183 // check for mid-edge points; reuse gv_verts 01184 gv_verts.resize(edges.size()); 01185 01186 // get the gventities 01187 result = MBI->tag_get_data(gvEntityHandle, edges, &gv_verts[0]); 01188 if (MB_SUCCESS != result) { 01189 std::cout << "Couldn't get gv edge data." << std::endl; 01190 return result; 01191 } 01192 for (rit = edges.begin(), i = 0; rit != edges.end(); rit++, i++) { 01193 assert(NULL != gv_verts[i]); 01194 int index = gv_verts[i]->get_index(dual_surf); 01195 assert(index >= 0); 01196 if (NULL != gv_verts[i]->gvizPoints[index] && 01197 gv_verts[i]->vtkEntityIds[index+2] == -1) { 01198 get_graphpoint_pos(gv_verts[i]->gvizPoints[index], dum_pos); 01199 if (useGraphviz) { 01200 sprintf(dum_str, "%d,%d", (int)(dum_pos[0]*x_xform), (int)(dum_pos[1]*y_xform)); 01201 agxset(gv_verts[i]->gvizPoints[index], asym_pos->index, dum_str); 01202 } 01203 01204 gv_verts[i]->vtkEntityIds[index+2] = 01205 points->InsertNextPoint(dum_pos[0], dum_pos[1], dum_pos[2]); 01206 vert_gv_map[*rit] = gv_verts[i]; 01207 } 01208 } 01209 01210 return MB_SUCCESS; 01211 } 01212 01213 ErrorCode DrawDual::get_xform(EntityHandle dual_surf, Agsym_t *asym_pos, 01214 double &x_xform, double &y_xform) 01215 { 01216 Range face_verts, face_verts_dum; 01217 01218 x_xform = y_xform = 1.0; 01219 return MB_SUCCESS; 01220 01221 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, NULL, NULL, 01222 &face_verts_dum, &face_verts, NULL); 01223 if (MB_SUCCESS != result) return result; 01224 01225 // get the gventities 01226 std::vector<GVEntity*> gv_verts; 01227 gv_verts.resize(face_verts.size()); 01228 result = MBI->tag_get_data(gvEntityHandle, face_verts, &gv_verts[0]); 01229 if (MB_SUCCESS != result) return result; 01230 01231 // find a vertex with non-zero x, y coordinates 01232 for (std::vector<GVEntity*>::iterator vit = gv_verts.begin(); vit != gv_verts.end(); vit++) { 01233 assert(NULL != *vit); 01234 int index = (*vit)->get_index(dual_surf); 01235 assert(index >= 0 && NULL != (*vit)->gvizPoints[index]); 01236 // get new point position 01237 double dum_pos[3]; 01238 get_graphpoint_pos((*vit)->gvizPoints[index], dum_pos); 01239 if (dum_pos[0] != 0 && dum_pos[1] != 0) { 01240 // get old point position, which is set on attribute 01241 char *agx = agxget((*vit)->gvizPoints[index], asym_pos->index); 01242 int x = -1, y = -1; 01243 sscanf(agx, "%d,%d", &x, &y); 01244 if (0 == x || 0 == y) continue; 01245 01246 // ok, non-zeros in all quantities, can compute xform now as old over new 01247 x_xform = ((double)x) / dum_pos[0]; 01248 y_xform = ((double)y) / dum_pos[1]; 01249 01250 return MB_SUCCESS; 01251 } 01252 } 01253 01254 // didn't get any usable data; set to unity 01255 x_xform = 1.0; 01256 y_xform = 1.0; 01257 01258 if (my_debug) std::cout << "Didn't find transform." << std::endl; 01259 01260 return MB_FAILURE; 01261 } 01262 01263 vtkPolyData *DrawDual::get_polydata(SheetDiagramPopup *this_sdpopup) 01264 { 01265 vtkRendererCollection *rcoll = this_sdpopup->sheet_diagram()->GetRenderWindow()->GetRenderers(); 01266 assert(rcoll != NULL); 01267 rcoll->InitTraversal(); 01268 vtkActorCollection *acoll = rcoll->GetNextItem()->GetActors(); 01269 assert(NULL != acoll); 01270 acoll->InitTraversal(); 01271 return vtkPolyData::SafeDownCast(acoll->GetNextActor()->GetMapper()->GetInput()); 01272 } 01273 01274 void DrawDual::get_clean_pd(EntityHandle dual_surf, 01275 SheetDiagramPopup *&this_sdpopup, 01276 vtkPolyData *&pd) 01277 { 01278 if (NULL != this_sdpopup) { 01279 ErrorCode result = reset_drawing_data(dual_surf); 01280 if (MB_SUCCESS != result) { 01281 std::cerr << "Trouble resetting drawing data for sheet." << std::endl; 01282 } 01283 } 01284 01285 if (NULL == this_sdpopup) { 01286 vtkRenderer *this_ren = vtkRenderer::New(); 01287 pd = vtkPolyData::New(); 01288 01289 const bool twod = false; 01290 01291 if (twod) { 01292 01293 // the easy route - just make new stuff 01294 vtkPolyDataMapper2D *this_mapper = vtkPolyDataMapper2D::New(); 01295 // vtkPolyDataMapper *this_mapper = vtkPolyDataMapper::New(); 01296 01297 // need to set a coordinate system for this window, so that display coordinates 01298 // are re-normalized to window size 01299 vtkCoordinate *this_coord = vtkCoordinate::New(); 01300 this_coord->SetCoordinateSystemToWorld(); 01301 this_mapper->SetTransformCoordinate(this_coord); 01302 this_coord->FastDelete(); 01303 01304 this_mapper->ScalarVisibilityOn(); 01305 this_mapper->SetLookupTable(vtkMOABUtils::lookupTable); 01306 this_mapper->UseLookupTableScalarRangeOn(); 01307 this_mapper->SetScalarModeToUseCellData(); 01308 01309 // put an edge extractor before the mapper 01310 vtkExtractEdges *ee = vtkExtractEdges::New(); 01311 ee->SetInput(pd); 01312 this_mapper->SetInput(ee->GetOutput()); 01313 this_mapper->SetInput(pd); 01314 ee->FastDelete(); 01315 01316 vtkActor2D *this_actor = vtkActor2D::New(); 01317 this_actor->PickableOn(); 01318 vtkMOABUtils::propSetMap[this_actor] = dual_surf; 01319 // vtkActor *this_actor = vtkActor::New(); 01320 this_actor->SetMapper(this_mapper); 01321 this_actor->GetProperty()->SetLineWidth(2.0); 01322 this_ren->AddActor(this_actor); 01323 this_actor->FastDelete(); 01324 this_mapper->FastDelete(); 01325 } 01326 else { 01327 // the easy route - just make new stuff 01328 vtkPolyDataMapper *this_mapper = vtkPolyDataMapper::New(); 01329 01330 // need to set a coordinate system for this window, so that display coordinates 01331 // are re-normalized to window size 01332 // vtkCoordinate *this_coord = vtkCoordinate::New(); 01333 // this_coord->SetCoordinateSystemToWorld(); 01334 // this_mapper->SetTransformCoordinate(this_coord); 01335 this_mapper->ScalarVisibilityOn(); 01336 //this_mapper->SetLookupTable(vtkMOABUtils::lookupTable); 01337 this_mapper->UseLookupTableScalarRangeOn(); 01338 this_mapper->SetScalarModeToUseCellData(); 01339 01340 // put an edge extractor before the mapper 01341 // vtkExtractEdges *ee = vtkExtractEdges::New(); 01342 // ee->SetInput(pd); 01343 // this_mapper->SetInput(ee->GetOutput()); 01344 this_mapper->SetInput(pd); 01345 01346 vtkActor *this_actor = vtkActor::New(); 01347 this_actor->PickableOn(); 01348 vtkMOABUtils::propSetMap[this_actor] = dual_surf; 01349 this_actor->SetMapper(this_mapper); 01350 this_mapper->FastDelete(); 01351 this_actor->GetProperty()->SetLineWidth(2.0); 01352 01353 double red, green, blue; 01354 int dum; 01355 vtkMOABUtils::get_colors(dual_surf, vtkMOABUtils::totalColors, dum, 01356 red, green, blue); 01357 vtkProperty *this_property = this_actor->GetProperty(); 01358 this_property->SetColor(red, green, blue); 01359 01360 this_ren->AddActor(this_actor); 01361 this_actor->FastDelete(); 01362 } 01363 01364 vtkCamera *camera = vtkCamera::New(); 01365 camera->SetPosition(CENT_X, CENT_Y, 3); 01366 camera->SetFocalPoint(CENT_X, CENT_Y, 0); 01367 camera->SetViewUp(0,1,0); 01368 01369 this_sdpopup = new SheetDiagramPopup(); 01370 if (my_debug) { 01371 //this_sdpopup->sheet_diagram()->GetRenderWindow()->DebugOn(); 01372 } 01373 01374 this_sdpopup->sheet_diagram()->GetRenderWindow()->AddRenderer(this_ren); 01375 this_sdpopup->sheet_diagram()->GetRenderWindow()->SetSize(SHEET_WINDOW_SIZE, SHEET_WINDOW_SIZE); 01376 01377 this_sdpopup->show(); 01378 01379 this_ren->ResetCamera(); 01380 this_ren->FastDelete(); 01381 01382 camera->Delete(); 01383 01384 return; 01385 } 01386 01387 // trace back until we find the dataset 01388 pd = get_polydata(this_sdpopup); 01389 assert(NULL != pd); 01390 // re-initialize the data, then we're done 01391 pd->Initialize(); 01392 } 01393 01394 ErrorCode DrawDual::construct_graphviz_data(EntityHandle dual_surf, 01395 Range &dcells, Range &dedges, 01396 Range &dverts, Range &face_verts, 01397 Range &loop_edges) 01398 { 01399 // gather/construct the data for graphviz 01400 01401 // get which drawring we're referring to 01402 char dum_str[80]; 01403 GraphWindows &this_gw = surfDrawrings[dual_surf]; 01404 if (useGraphviz && NULL == this_gw.gvizGraph) { 01405 // allocate a new graph and create a new window 01406 sprintf(dum_str, "%d", (int)MBI->id_from_handle(dual_surf)); 01407 this_gw.gvizGraph = agopen(dum_str, AGRAPH); 01408 if (this_gw.gvizGraph == NULL) return MB_FAILURE; 01409 } 01410 01411 // for each vertex, allocate a graphviz point if it doesn't already have one 01412 Agsym_t *asym_pos = (useGraphviz ? get_asym(dual_surf, 0, "pos") : NULL); 01413 01414 ErrorCode result = construct_graphviz_points(dual_surf, dverts, asym_pos); RR; 01415 01416 // for each edge, allocate a graphviz edge if it doesn't already have one 01417 result = construct_graphviz_edges(dual_surf, dedges, face_verts, asym_pos); RR; 01418 01419 return result; 01420 } 01421 01422 ErrorCode DrawDual::construct_graphviz_edges(EntityHandle dual_surf, 01423 Range &dedges, 01424 Range &loop_verts, 01425 Agsym_t *asym_pos) 01426 { 01427 const EntityHandle *connect; 01428 int num_connect; 01429 ErrorCode result = MB_SUCCESS; 01430 char dum_str[80]; 01431 GraphWindows &this_gw = surfDrawrings[dual_surf]; 01432 //Agsym_t *asym_weight = (useGraphviz ? get_asym(dual_surf, 1, "weight") : NULL), 01433 //*asym_len = (useGraphviz ? get_asym(dual_surf, 1, "len") : NULL); 01434 std::vector<GVEntity*> dedge_gv(dedges.size()); 01435 result = MBI->tag_get_data(gvEntityHandle, dedges, &dedge_gv[0]); RR; 01436 GVEntity *dvert_gv[2]; 01437 01438 Range::iterator rit; 01439 int i; 01440 for (rit = dedges.begin(), i = 0; rit != dedges.end(); rit++, i++) { 01441 01442 // get the DEdgeGVEdge for this dual edge 01443 GVEntity *this_gv = dedge_gv[i]; 01444 if (NULL == this_gv) { 01445 this_gv = new GVEntity(); 01446 this_gv->moabEntity = *rit; 01447 result = MBI->tag_set_data(gvEntityHandle, &(*rit), 1, &this_gv); 01448 if (MB_SUCCESS != result) return result; 01449 dedge_gv[i] = this_gv; 01450 } 01451 01452 // get the index for this dual surface 01453 int dsindex = this_gv->get_index(dual_surf); 01454 assert(dsindex > -10); 01455 if (dsindex < 0) { 01456 // need to make a graphviz entity; get the graphviz nodes, then make the edge 01457 // next sheet index is negative of value returned from get_index 01458 dsindex = -dsindex - 1; 01459 this_gv->dualSurfs[dsindex] = dual_surf; 01460 01461 result = MBI->get_connectivity(*rit, connect, num_connect); 01462 if (MB_SUCCESS != result) return result; 01463 result = MBI->tag_get_data(gvEntityHandle, connect, 2, dvert_gv); 01464 if (MB_SUCCESS != result) return result; 01465 assert(NULL != dvert_gv[0] && NULL != dvert_gv[1]); 01466 int index0 = dvert_gv[0]->get_index(dual_surf); 01467 int index1 = dvert_gv[1]->get_index(dual_surf); 01468 assert(index0 >= 0 && index1 >= 0); 01469 sprintf(dum_str, "%d", (int)*rit); 01470 01471 // first, check to see if it's degenerate; if so, add a mid-pt 01472 Range tmp_edges; 01473 result = MBI->get_adjacencies(connect, 2, 1, false, tmp_edges); 01474 if (MB_SUCCESS == result && tmp_edges.size() > 1) { 01475 if (useGraphviz) { 01476 // add a graphviz pt for the edge 01477 Agnode_t *mid_gvpt = agnode(this_gw.gvizGraph, dum_str); 01478 sprintf(dum_str, "%up", (int)MBI->id_from_handle(*rit)); 01479 this_gv->gvizPoints[dsindex] = mid_gvpt; 01480 01481 Agedge_t *this_gved = agedge(this_gw.gvizGraph, 01482 (Agnode_t*)dvert_gv[0]->gvizPoints[index0], 01483 (Agnode_t*)mid_gvpt); 01484 this_gv->gvizEdges[dsindex] = this_gved; 01485 01486 this_gved = agedge(this_gw.gvizGraph, 01487 mid_gvpt, 01488 (Agnode_t*)dvert_gv[1]->gvizPoints[index1]); 01489 this_gv->gvizEdges[dsindex+2] = this_gved; 01490 } 01491 else { 01492 // add a vertex for the edge 01493 sprintf(dum_str, "%up", (unsigned int) MBI->id_from_handle(*rit)); 01494 EntityHandle mid_vert, edge1, edge2, edge_verts[2]; 01495 double dum_pos[] = {CENT_X, CENT_Y, 0.0}; 01496 result = MBI->create_vertex(dum_pos, mid_vert); 01497 if (MB_SUCCESS != result) return result; 01498 this_gv->gvizPoints[dsindex] = (void*)mid_vert; 01499 01500 edge_verts[0] = (EntityHandle) dvert_gv[0]->gvizPoints[index0]; 01501 edge_verts[1] = (EntityHandle) mid_vert; 01502 result = MBI->create_element(MBEDGE, edge_verts, 2, edge1); 01503 if (MB_SUCCESS != result) return result; 01504 this_gv->gvizEdges[dsindex] = (void*)edge1; 01505 result = MBI->tag_set_data(dualCurveTagHandle, &edge1, 1, &(*rit)); 01506 if (MB_SUCCESS != result) return result; 01507 01508 edge_verts[0] = (EntityHandle) mid_vert; 01509 edge_verts[1] = (EntityHandle) dvert_gv[1]->gvizPoints[index1]; 01510 result = MBI->create_element(MBEDGE, edge_verts, 2, edge2); 01511 if (MB_SUCCESS != result) return result; 01512 this_gv->gvizEdges[dsindex+2] = (void*)edge2; 01513 result = MBI->tag_set_data(dualCurveTagHandle, &edge2, 1, &(*rit)); 01514 if (MB_SUCCESS != result) return result; 01515 } 01516 } 01517 else { 01518 if (useGraphviz) { 01519 Agedge_t *this_gved = agedge(this_gw.gvizGraph, 01520 (Agnode_t*)dvert_gv[0]->gvizPoints[index0], 01521 (Agnode_t*)dvert_gv[1]->gvizPoints[index1]); 01522 this_gv->gvizEdges[dsindex] = this_gved; 01523 } 01524 else { 01525 EntityHandle edge_verts[2], edge1; 01526 edge_verts[0] = (EntityHandle) dvert_gv[0]->gvizPoints[index0]; 01527 edge_verts[1] = (EntityHandle) dvert_gv[1]->gvizPoints[index1]; 01528 result = MBI->create_element(MBEDGE, edge_verts, 2, edge1); 01529 if (MB_SUCCESS != result) return result; 01530 this_gv->gvizEdges[dsindex] = (void*) edge1; 01531 result = MBI->tag_set_data(dualCurveTagHandle, &edge1, 1, &(*rit)); 01532 if (MB_SUCCESS != result) return result; 01533 } 01534 } 01535 01536 // check to see if it's an interior edge connected to the loop, and 01537 // adjust its weight and length if so 01538 Range::iterator firstv = loop_verts.find(connect[0]), 01539 lastv = loop_verts.find(connect[1]), 01540 endloop = loop_verts.end(); 01541 if ((firstv == endloop && lastv != endloop) || 01542 (firstv != endloop && lastv == endloop)) { 01543 // agxset(this_gved, asym_weight->index, "10"); 01544 // agxset(this_gved, asym_len->index, "0.1"); 01545 } 01546 } 01547 01548 // note to self: check for self-intersections here, to handle 2 instances 01549 // of edge in drawring 01550 } 01551 01552 return result; 01553 } 01554 01555 ErrorCode DrawDual::construct_graphviz_points(EntityHandle dual_surf, 01556 Range &dverts, 01557 Agsym_t *asym_pos) 01558 { 01559 Range::iterator rit; 01560 int i; 01561 ErrorCode result = MB_SUCCESS; 01562 GraphWindows &this_gw = surfDrawrings[dual_surf]; 01563 char dum_str[80]; 01564 std::vector<GVEntity*> dvert_gv(dverts.size()); 01565 result = MBI->tag_get_data(gvEntityHandle, dverts, &dvert_gv[0]); 01566 01567 for (rit = dverts.begin(), i = 0; rit != dverts.end(); rit++, i++) { 01568 01569 // get the DVertexPoint for this dual vertex 01570 GVEntity *this_gv = dvert_gv[i]; 01571 if (NULL == this_gv) { 01572 this_gv = new GVEntity(); 01573 this_gv->moabEntity = *rit; 01574 result = MBI->tag_set_data(gvEntityHandle, &(*rit), 1, &this_gv); 01575 if (MB_SUCCESS != result) return result; 01576 dvert_gv[i] = this_gv; 01577 } 01578 01579 // get the index for this dual surface 01580 int dsindex = this_gv->get_index(dual_surf); 01581 assert(dsindex > -10); 01582 if (dsindex < 0) { 01583 dsindex = -dsindex - 1; 01584 if (useGraphviz) { 01585 // need to make a graphviz point 01586 sprintf(dum_str, "%u", (unsigned int)*rit); 01587 Agnode_t *this_gvpt = agnode(this_gw.gvizGraph, dum_str); 01588 if (NULL == this_gvpt) return MB_FAILURE; 01589 this_gv->gvizPoints[dsindex] = this_gvpt; 01590 } 01591 else { 01592 // use an MBVertex instead 01593 EntityHandle new_vertex; 01594 double dum_pos[] = {CENT_X, CENT_Y, 0.0}; 01595 result = MBI->create_vertex(dum_pos, new_vertex); 01596 if (MB_SUCCESS != result) return result; 01597 this_gv->gvizPoints[dsindex] = (void*) new_vertex; 01598 } 01599 01600 this_gv->dualSurfs[dsindex] = dual_surf; 01601 } 01602 } 01603 01604 return result; 01605 } 01606 01607 Agsym_t *DrawDual::get_asym(EntityHandle dual_surf, const int dim, 01608 const char *name, const char *def_val) 01609 { 01610 Agsym_t *asym = NULL; 01611 01612 if (0 == dim) { 01613 asym = agfindattr(surfDrawrings[dual_surf].gvizGraph->proto->n, 01614 const_cast<char*>(name)); 01615 if (NULL == asym) { 01616 if (NULL == def_val) asym = agnodeattr(surfDrawrings[dual_surf].gvizGraph, const_cast<char*>(name), (char*)""); 01617 01618 else asym = agnodeattr(surfDrawrings[dual_surf].gvizGraph, const_cast<char*>(name), 01619 const_cast<char*>(def_val)); 01620 } 01621 } 01622 else { 01623 asym = agfindattr(surfDrawrings[dual_surf].gvizGraph->proto->e, 01624 const_cast<char*>(name)); 01625 if (NULL == asym) { 01626 if (NULL == def_val) asym = agedgeattr(surfDrawrings[dual_surf].gvizGraph, const_cast<char*>(name), (char*)""); 01627 01628 else asym = agedgeattr(surfDrawrings[dual_surf].gvizGraph, const_cast<char*>(name), 01629 const_cast<char*>(def_val)); 01630 } 01631 } 01632 01633 assert(NULL != asym); 01634 return asym; 01635 } 01636 01637 ErrorCode DrawDual::compute_fixed_points(EntityHandle dual_surf, Range &dverts, 01638 Range &face_verts, Range &loop_edges) 01639 { 01640 std::vector<std::vector<EntityHandle> > loops; 01641 Range new_face_verts(face_verts); 01642 01643 while (!new_face_verts.empty()) { 01644 // get the next first vertex on the loop 01645 EntityHandle this_v = *new_face_verts.begin(); 01646 EntityHandle first_v = 0, last_v = 0; 01647 std::vector<EntityHandle> loop_vs; 01648 Range temp_face_verts; 01649 01650 if (new_face_verts.size() == 2) { 01651 // quick way to do this, assuming both vertices are on the loop 01652 loop_vs.push_back(*new_face_verts.begin()); 01653 loop_vs.push_back(*new_face_verts.rbegin()); 01654 temp_face_verts.insert(*new_face_verts.begin()); 01655 temp_face_verts.insert(*new_face_verts.rbegin()); 01656 this_v = first_v; 01657 } 01658 01659 while (this_v != first_v) { 01660 if (0 == first_v) first_v = this_v; 01661 01662 // put this vertex on the loop, then get the next one 01663 loop_vs.push_back(this_v); 01664 temp_face_verts.insert(this_v); 01665 EntityHandle temp_v = this_v; 01666 this_v = vtkMOABUtils::dualTool->next_loop_vertex(last_v, this_v, dual_surf); 01667 assert(0 != this_v); 01668 last_v = temp_v; 01669 } 01670 01671 // save this vector in the map 01672 loops.push_back(loop_vs); 01673 01674 // ok, we've got them all; first, remove them from face_verts 01675 Range temp_range = subtract( new_face_verts, temp_face_verts); 01676 new_face_verts.swap(temp_range); 01677 } 01678 01679 // now compute vertex coordinates for each loop 01680 Agsym_t *asym_pos, *asym_pin; 01681 if (useGraphviz) { 01682 asym_pos = get_asym(dual_surf, 0, "pos"); 01683 asym_pin = get_asym(dual_surf, 0, "pin", "false"); 01684 } 01685 01686 if (loops.empty()) 01687 return compute_pillow_fixed_points(dual_surf, face_verts, loop_edges); 01688 01689 char tmp_pos[80]; 01690 int loop_num, num_loops = loops.size(); 01691 std::vector<std::vector<EntityHandle> >::iterator mit; 01692 if (my_debug) std::cout << "Loop points: " << std::endl; 01693 for (mit = loops.begin(), loop_num = 0; mit != loops.end(); mit++, loop_num++) { 01694 01695 // if we have more than two loops, let graphviz compute the best position 01696 if (num_loops > 2 && loop_num > 0) break; 01697 01698 // now, go around the loop, assigning vertex positions; assume a 6-inch diameter circle 01699 // (at 72 pts/inch) 01700 unsigned int loop_size = (*mit).size(); 01701 double angle = (2.0*acos(-1.0))/((double)loop_size); 01702 int xpos_pts, ypos_pts; 01703 01704 for (unsigned int i = 0; i < loop_size; i++) { 01705 // get the position 01706 get_loop_vertex_pos(i, loop_num, num_loops, angle, xpos_pts, ypos_pts); 01707 01708 // now set that position on the node 01709 GVEntity *this_gv; 01710 ErrorCode result = MBI->tag_get_data(gvEntityHandle, &((*mit)[i]), 1, &this_gv); RR; 01711 01712 int index = this_gv->get_index(dual_surf); 01713 assert(index >= 0); 01714 01715 if (useGraphviz) { 01716 Agnode_t *this_gpt = (Agnode_t*)this_gv->gvizPoints[index]; 01717 01718 // set position and pin attributes for the node 01719 sprintf(tmp_pos, "%d,%d!", xpos_pts, ypos_pts); 01720 agxset(this_gpt, asym_pos->index, tmp_pos); 01721 agxset(this_gpt, asym_pin->index, (char*)"true"); 01722 01723 // also try setting them in the data structure directly 01724 ND_coord_i(this_gpt).x = xpos_pts; 01725 ND_coord_i(this_gpt).y = ypos_pts; 01726 ND_pinned(this_gpt) = true; 01727 } 01728 else { 01729 EntityHandle this_vert = (EntityHandle) this_gv->gvizPoints[index]; 01730 double dum_pos[] = {xpos_pts, ypos_pts, 0.0}; 01731 result = MBI->set_coords(&this_vert, 1, dum_pos); 01732 if (MB_SUCCESS != result) return result; 01733 } 01734 01735 if (my_debug) std::cout << "Point " << MBI->id_from_handle((*mit)[i]) 01736 << ": x = " << xpos_pts << ", y = " << ypos_pts << std::endl; 01737 } 01738 01739 if (loop_size == 2) { 01740 // 2-vertex loop: place mid-edge nodes too 01741 assert(loop_edges.size() == 2); 01742 01743 int offset = 0; 01744 // put one at pi/2, the other at 3pi/2 01745 for (Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); rit++) { 01746 GVEntity *this_gv; 01747 ErrorCode result = MBI->tag_get_data(gvEntityHandle, &(*rit), 1, &this_gv); 01748 if (MB_SUCCESS != result) continue; 01749 int index = this_gv->get_index(dual_surf); 01750 assert(index >= 0 && this_gv->gvizPoints[index] != NULL); 01751 get_loop_vertex_pos(1+2*offset, loop_num, num_loops, .5*angle, xpos_pts, ypos_pts); 01752 01753 if (useGraphviz) { 01754 Agnode_t *this_gpt = (Agnode_t*)this_gv->gvizPoints[index]; 01755 sprintf(tmp_pos, "%d,%d!", xpos_pts, ypos_pts); 01756 agxset(this_gpt, asym_pos->index, tmp_pos); 01757 agxset(this_gpt, asym_pin->index, (char*)"true"); 01758 01759 // also try setting them in the data structure directly 01760 ND_coord_i(this_gpt).x = xpos_pts; 01761 ND_coord_i(this_gpt).y = ypos_pts; 01762 ND_pinned(this_gpt) = true; 01763 } 01764 else { 01765 EntityHandle this_vert = (EntityHandle) this_gv->gvizPoints[index]; 01766 double dum_pos[] = {xpos_pts, ypos_pts, 0.0}; 01767 result = MBI->set_coords(&this_vert, 1, dum_pos); 01768 if (MB_SUCCESS != result) return result; 01769 } 01770 01771 if (my_debug) std::cout << "Edge point for edge " << MBI->id_from_handle(*rit) 01772 << ": x = " << xpos_pts << ", y = " << ypos_pts << std::endl; 01773 offset += 1; 01774 } 01775 } 01776 } 01777 01778 return MB_SUCCESS; 01779 } 01780 01781 ErrorCode DrawDual::compute_pillow_fixed_points(EntityHandle dual_surf, 01782 Range &face_verts, 01783 Range &face_edges) 01784 { 01785 // find the points we'll call fixed for this sheet 01786 // first, get the chords, in an ordered list 01787 std::vector<EntityHandle> chords; 01788 ErrorCode result = MBI->get_child_meshsets(dual_surf, chords); 01789 if (MB_SUCCESS != result) return result; 01790 01791 // if there are only two, don't bother, they'll get fixed up later 01792 if (chords.size() <= 2) return MB_SUCCESS; 01793 01794 // if more, get the first two vertices on the first chord, then a 3rd 01795 // from another that's not in the 1st set 01796 Range tmp_range; 01797 result = MBI->get_entities_by_dimension(chords[0], 1, tmp_range); RR; 01798 const EntityHandle *connect; 01799 int num_connect; 01800 result = MBI->get_connectivity(*tmp_range.begin(), connect, num_connect); RR; 01801 face_verts.insert(connect[0]); 01802 face_verts.insert(connect[1]); 01803 tmp_range.clear(); 01804 Range tmp_verts; 01805 result = MBI->get_entities_by_dimension(chords[1], 1, tmp_range); RR; 01806 result = MBI->get_adjacencies(tmp_range, 0, false, tmp_verts, Interface::UNION); RR; 01807 tmp_verts = subtract( tmp_verts, face_verts); 01808 assert(!tmp_verts.empty()); 01809 face_verts.insert(*tmp_verts.begin()); 01810 01811 // should have 3; put at equal angles round the circle 01812 // double PI = acos(-1.0); 01813 double dum_pos[9]; 01814 dum_pos[0] = CENT_X; 01815 dum_pos[1] = RAD_PTS; 01816 dum_pos[3] = CENT_X - RAD_PTS*cos(PI/6.0); 01817 dum_pos[4] = CENT_Y - RAD_PTS*sin(PI/6.0); 01818 dum_pos[6] = CENT_X + RAD_PTS*cos(PI/6.0); 01819 dum_pos[7] = dum_pos[4]; 01820 dum_pos[2] = dum_pos[5] = dum_pos[8] = 0.0; 01821 01822 std::vector<EntityHandle> graph_points(face_verts.size()); 01823 get_graph_points(face_verts, dual_surf, (void**) &graph_points[0]); 01824 result = MBI->set_coords(&graph_points[0], face_verts.size(), dum_pos); RR; 01825 01826 return result; 01827 } 01828 01829 void DrawDual::get_loop_vertex_pos(unsigned int vert_num, 01830 unsigned int loop_num, 01831 unsigned int num_loops, 01832 double angle, int &xpos_pts, int &ypos_pts) 01833 { 01834 double this_angle = vert_num * angle; 01835 if (num_loops > 2 && loop_num > 0) { 01836 xpos_pts = 0; 01837 ypos_pts = 0; 01838 return; 01839 } 01840 01841 xpos_pts = (int) (((double)RAD_PTS) * cos(this_angle)); 01842 ypos_pts = (int) (((double)RAD_PTS) * sin(this_angle)); 01843 01844 if (loop_num > 0) { 01845 xpos_pts /= 6; 01846 ypos_pts /= 6; 01847 } 01848 01849 xpos_pts += CENT_X; 01850 ypos_pts += CENT_Y; 01851 } 01852 01853 ErrorCode DrawDual::draw_labels(EntityHandle dual_surf, vtkPolyData *pd, 01854 vtkPolyData *new_pd) 01855 { 01856 // get the renderer you'll be writing the text to 01857 GraphWindows &this_gw = surfDrawrings[dual_surf]; 01858 vtkRenderer *ren = 01859 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->GetRenderers()->GetFirstRenderer(); 01860 01861 // sheet id first 01862 char set_name[CATEGORY_TAG_SIZE]; 01863 int dum; 01864 ErrorCode result = vtkMOABUtils::MBI->tag_get_data(vtkMOABUtils::globalId_tag(), 01865 &dual_surf, 1, &dum); 01866 if (MB_SUCCESS != result) return result; 01867 sprintf(set_name, "%d\n", dum); 01868 01869 // create a text actor 01870 vtkTextActor *text_actor = vtkTextActor::New(); 01871 vtkMOABUtils::propSetMap[text_actor] = dual_surf; 01872 text_actor->SetInput(set_name); 01873 01874 // compute proper position for string 01875 double LABEL_FRACTION = .90; 01876 vtkCoordinate *this_pos = text_actor->GetPositionCoordinate(); 01877 this_pos->SetCoordinateSystemToWorld(); 01878 this_pos->SetValue(LABEL_FRACTION*RAD_PTS+CENT_X, LABEL_FRACTION*RAD_PTS+CENT_Y, 0); 01879 text_actor->GetTextProperty()->SetColor(1.0, 1.0, 1.0); 01880 text_actor->GetTextProperty()->BoldOn(); 01881 ren->AddActor(text_actor); 01882 text_actor->FastDelete(); 01883 01884 // now vertex ids 01885 01886 // first, get new polydata's with interior and loop points 01887 label_interior_verts(dual_surf, pd, ren); 01888 01889 // create other sheet labels and add to renderer 01890 vtkLabeledDataMapper *os_labels = vtkLabeledDataMapper::New(); 01891 os_labels->SetInput(new_pd); 01892 os_labels->SetLabelModeToLabelFieldData(); 01893 os_labels->SetLabelFormat("___/%g"); 01894 vtkActor2D *lda2 = vtkActor2D::New(); 01895 lda2->SetMapper(os_labels); 01896 ren->AddActor(lda2); 01897 01898 // set label actor to be non-pickable 01899 lda2->PickableOff(); 01900 01901 os_labels->FastDelete(); 01902 lda2->FastDelete(); 01903 01904 return MB_SUCCESS; 01905 } 01906 01907 ErrorCode DrawDual::label_other_sheets(EntityHandle dual_surf, 01908 vtkPolyData *pd, 01909 vtkPolyData *&new_pd) 01910 { 01911 // gather the chord sets 01912 std::vector<EntityHandle> chords, chord_surfs; 01913 Range dedges, dverts, dverts_loop; 01914 ErrorCode result = MBI->get_child_meshsets(dual_surf, chords); 01915 if (MB_SUCCESS != result || 0 == chords.size()) return result; 01916 01917 // start a new polydata with points for labeling 01918 new_pd = vtkPolyData::New(); 01919 vtkPoints *new_points = new_pd->GetPoints(); 01920 if (NULL == new_points) { 01921 new_points = vtkPoints::New(); 01922 new_pd->SetPoints(new_points); 01923 new_points->Delete(); 01924 } 01925 new_pd->Allocate(); 01926 vtkIntArray *id_array = vtkIntArray::New(); 01927 id_array->SetName("LoopVertexIds"); 01928 01929 // for each chord: 01930 ErrorCode tmp_result; 01931 std::vector<EntityHandle>::iterator vit1; 01932 Range::iterator rit; 01933 std::vector<GVEntity*> gv_edges; 01934 for (vit1 = chords.begin(); vit1 != chords.end(); vit1++) { 01935 01936 // get a color for this chord; make it the color of the "other" sheet, unless 01937 // it's this sheet, in which case it's black 01938 EntityHandle color_set = other_sheet(*vit1, dual_surf); 01939 double red, green, blue; 01940 int global_id; 01941 if (0 == color_set) red = green = blue = 0.0; 01942 else 01943 vtkMOABUtils::get_colors(color_set, vtkMOABUtils::totalColors, global_id, 01944 red, green, blue); 01945 01946 // create a series of edges in the original pd 01947 dedges.clear(); dverts.clear(); dverts_loop.clear(); 01948 tmp_result = vtkMOABUtils::dualTool->get_dual_entities(*vit1, NULL, &dedges, &dverts, &dverts_loop, 01949 NULL); 01950 if (MB_SUCCESS != tmp_result) { 01951 result = tmp_result; 01952 continue; 01953 } 01954 01955 gv_edges.reserve(dedges.size()); 01956 tmp_result = MBI->tag_get_data(gvEntityHandle, dedges, &gv_edges[0]); 01957 if (MB_SUCCESS != tmp_result) { 01958 result = tmp_result; 01959 continue; 01960 } 01961 01962 int edge_vtk_vs[2]; 01963 GVEntity *gv_verts[2]; 01964 const EntityHandle *edge_vs; 01965 int num_edge_vs; 01966 int edge_num; 01967 int index; 01968 for (rit = dedges.begin(), edge_num = 0; rit != dedges.end(); rit++, edge_num++) { 01969 tmp_result = MBI->get_connectivity(*rit, edge_vs, num_edge_vs); 01970 if (MB_SUCCESS != tmp_result) { 01971 result = tmp_result; 01972 continue; 01973 } 01974 // get the gventities 01975 tmp_result = MBI->tag_get_data(gvEntityHandle, edge_vs, 2, gv_verts); 01976 if (MB_SUCCESS != tmp_result) { 01977 result = tmp_result; 01978 continue; 01979 } 01980 for (int i = 0; i < 2; i++) { 01981 index = gv_verts[i]->get_index(dual_surf); 01982 assert(index >= 0 && NULL != gv_verts[i]->gvizPoints[index]); 01983 edge_vtk_vs[i] = gv_verts[i]->vtkEntityIds[index]; 01984 01985 // look for loop points, and add label if there are any 01986 if (dverts_loop.find(edge_vs[i]) != dverts_loop.end()) { 01987 double dum_pos[3]; 01988 get_graphpoint_pos(gv_verts[i]->gvizPoints[index], dum_pos); 01989 new_points->InsertNextPoint(dum_pos[0], dum_pos[1], dum_pos[2]); 01990 id_array->InsertNextValue(global_id); 01991 } 01992 } 01993 } 01994 } 01995 01996 // assign id data to the pd for drawing other sheet labels 01997 new_pd->GetPointData()->AddArray(id_array); 01998 new_pd->GetPointData()->SetScalars(id_array); 01999 new_pd->GetPointData()->SetActiveAttribute("LoopVertexIds", 0); 02000 id_array->FastDelete(); 02001 02002 return result; 02003 } 02004 02005 void DrawDual::label_interior_verts(EntityHandle dual_surf, 02006 vtkPolyData *pd, 02007 vtkRenderer *ren) 02008 { 02009 // this function originally designed to filter out interior vertices, 02010 // extract them and add data just to them; however, there isn't a filter 02011 // in vtk to extract just vertices, so we label all the vertices here 02012 // 02013 // get the cells and vertices on this dual surface 02014 Range dcells, dverts, face_verts; 02015 02016 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, &dcells, NULL, 02017 &dverts, &face_verts, NULL); 02018 if (MB_SUCCESS != result) return; 02019 02020 // get the GVentity for vertices, so we can find the vtk point on the sheet drawing 02021 std::vector<GVEntity*> gv_ents; 02022 gv_ents.reserve(dverts.size()); 02023 result = MBI->tag_get_data(gvEntityHandle, dverts, &gv_ents[0]); 02024 if (MB_SUCCESS != result) { 02025 std::cout << "Failed to get GV entities." << std::endl; 02026 return; 02027 } 02028 02029 // get the ids of the primal entities of these vertices 02030 std::vector<int> vert_ids; 02031 result = get_primal_ids(dverts, vert_ids); 02032 if (MB_SUCCESS != result) { 02033 std::cout << "Failed to get primal ids." << std::endl; 02034 return; 02035 } 02036 02037 // now put those ids in the id list, and the vtk point ids 02038 // in an extract cells list 02039 vtkPolyData *label_pd = vtkPolyData::New(); 02040 vtkPoints *new_points = label_pd->GetPoints(); 02041 if (NULL == new_points) { 02042 new_points = vtkPoints::New(); 02043 label_pd->SetPoints(new_points); 02044 new_points->Delete(); 02045 } 02046 label_pd->Allocate(); 02047 vtkIntArray *int_ids = vtkIntArray::New(); 02048 int_ids->SetName("VertexIds"); 02049 int_ids->SetNumberOfValues(dverts.size()); 02050 02051 for (unsigned int i = 0; i != dverts.size(); i++) { 02052 int index = gv_ents[i]->get_index(dual_surf); 02053 assert(index >= 0); 02054 02055 // insert the primal entity id into a list for labeling 02056 double dum_pos[3]; 02057 get_graphpoint_pos(gv_ents[i]->gvizPoints[index], dum_pos); 02058 new_points->InsertNextPoint(dum_pos[0], dum_pos[1], dum_pos[2]); 02059 int_ids->InsertValue(i, vert_ids[i]); 02060 } 02061 02062 // make a new pd and copy the points from the last one 02063 label_pd->GetPointData()->AddArray(int_ids); 02064 label_pd->GetPointData()->SetActiveAttribute("VertexIds", 0); 02065 02066 vtkLabeledDataMapper *ldm = vtkLabeledDataMapper::New(); 02067 02068 ldm->SetInput(label_pd); 02069 ldm->SetLabelModeToLabelScalars(); 02070 ldm->SetLabelFormat("%g"); 02071 vtkActor2D *lda = vtkActor2D::New(); 02072 lda->SetMapper(ldm); 02073 ren->AddActor(lda); 02074 02075 // set label actor to be non-pickable 02076 lda->PickableOff(); 02077 02078 label_pd->FastDelete(); 02079 int_ids->FastDelete(); 02080 ldm->FastDelete(); 02081 lda->FastDelete(); 02082 } 02083 02084 EntityHandle DrawDual::other_sheet(const EntityHandle this_chord, 02085 const EntityHandle dual_surf) 02086 { 02087 static std::vector<EntityHandle> chord_surfs; 02088 EntityHandle val = dual_surf; 02089 ErrorCode result = MBI->get_parent_meshsets(this_chord, chord_surfs); 02090 if (MB_SUCCESS == result && !chord_surfs.empty()) { 02091 if (chord_surfs[0] != dual_surf) val = chord_surfs[0]; 02092 else if (chord_surfs.size() > 1 && chord_surfs[1] != dual_surf) 02093 val = chord_surfs[1]; 02094 else if (chord_surfs[0] == dual_surf && 02095 (chord_surfs.size() == 1 || chord_surfs[1] == dual_surf)) 02096 val = 0; 02097 } 02098 chord_surfs.clear(); 02099 return val; 02100 } 02101 02102 ErrorCode DrawDual::get_primal_ids(const Range &ents, std::vector<int> &ids) 02103 { 02104 // get the ids of these verts, equal to entity ids of their primal entities 02105 static std::vector<EntityHandle> primals; 02106 primals.resize(ents.size()); 02107 ids.resize(ents.size()); 02108 ErrorCode result = MBI->tag_get_data(dualEntityTagHandle, ents, &primals[0]); 02109 if (MB_SUCCESS != result) { 02110 for (unsigned int i = 0; i < ents.size(); i++) ids[i] = 0; 02111 } 02112 02113 result = MBI->tag_get_data(vtkMOABUtils::globalId_tag(), &primals[0], ents.size(), 02114 &ids[0]); 02115 if (MB_SUCCESS != result && MB_TAG_NOT_FOUND != result) 02116 std::cerr << "tag_get_data returned non-zero result in get_primal_ids." << std::endl; 02117 02118 for (unsigned int i = 0; i < ents.size(); i++) { 02119 if (0 == ids[i]) ids[i] = MBI->id_from_handle(primals[i]); 02120 } 02121 02122 return MB_SUCCESS; 02123 } 02124 02125 ErrorCode DrawDual::reset_drawn_sheets(Range *drawn_sheets) 02126 { 02127 ErrorCode result = MB_SUCCESS, tmp_result; 02128 for (std::map<EntityHandle,GraphWindows>::iterator mit = surfDrawrings.begin(); 02129 mit != surfDrawrings.end(); mit++) { 02130 if (NULL != drawn_sheets && 02131 NULL != (*mit).second.sheetDiagram) 02132 drawn_sheets->insert((*mit).first); 02133 tmp_result = (*mit).second.reset((*mit).first); 02134 if (MB_SUCCESS != tmp_result) result = tmp_result; 02135 } 02136 02137 return tmp_result; 02138 } 02139 02140 ErrorCode DrawDual::GraphWindows::reset(EntityHandle dual_surf) 02141 { 02142 ErrorCode result = MB_SUCCESS; 02143 02144 if (sheetDiagram) { 02145 result = gDrawDual->reset_drawing_data(dual_surf); 02146 delete sheetDiagram; 02147 sheetDiagram = NULL; 02148 } 02149 02150 if (pickActor) { 02151 pickActor->Delete(); 02152 pickActor = NULL; 02153 } 02154 02155 return result; 02156 } 02157 02158 ErrorCode DrawDual::reset_drawing_data(EntityHandle dual_surf) 02159 { 02160 // deleting a sheet drawing; reset the data on moab tags so it'll draw right 02161 // next time 02162 02163 // get the widget 02164 GraphWindows &this_gw = surfDrawrings[dual_surf]; 02165 02166 vtkRenderer *ren = 02167 this_gw.sheetDiagram->sheet_diagram()->GetRenderWindow()->GetRenderers()->GetFirstRenderer(); 02168 02169 vtkActorCollection *acoll = ren->GetActors(); 02170 assert(NULL != acoll); 02171 acoll->InitTraversal(); 02172 vtkActor *tmp_actor; 02173 Range saved_sets; 02174 02175 // get all actors, check sets in propSetMap; save set, remove actor from map 02176 while ((tmp_actor = acoll->GetNextItem())) { 02177 std::map<vtkProp*,EntityHandle>::iterator mit = 02178 vtkMOABUtils::propSetMap.find(tmp_actor); 02179 if (mit == vtkMOABUtils::propSetMap.end()) continue; 02180 saved_sets.insert((*mit).second); 02181 vtkMOABUtils::propSetMap.erase(mit); 02182 } 02183 02184 // for dual surface set: 02185 // get 0-, 1-, 2-cells, GVEntity for each 02186 Range tcells, all_cells; 02187 ErrorCode result = MBI->get_entities_by_type(dual_surf, MBPOLYGON, 02188 tcells); 02189 if (MB_SUCCESS != result) return result; 02190 result = MBI->get_adjacencies(tcells, 0, false, all_cells, Interface::UNION); 02191 if (MB_SUCCESS != result) return result; 02192 result = MBI->get_adjacencies(tcells, 1, false, all_cells, Interface::UNION); 02193 if (MB_SUCCESS != result) return result; 02194 02195 for (Range::const_reverse_iterator rit = all_cells.rbegin(); rit != all_cells.rend(); rit++) { 02196 // get the GVEntity 02197 GVEntity *gv_ent; 02198 result = MBI->tag_get_data(gvEntityHandle, &(*rit), 1, &gv_ent); 02199 if (MB_TAG_NOT_FOUND == result || 0 == gv_ent) continue; 02200 02201 // reset the data on this gv_ent for this dual surf 02202 int index = gv_ent->get_index(dual_surf); 02203 if (index >= 0) gv_ent->reset(index); 02204 } 02205 02206 if (this_gw.gvizGraph) { 02207 free(this_gw.gvizGraph); 02208 this_gw.gvizGraph = NULL; 02209 } 02210 02211 if (this_gw.pickActor) { 02212 this_gw.pickActor->Delete(); 02213 this_gw.pickActor = NULL; 02214 } 02215 02216 return MB_SUCCESS; 02217 } 02218 02219 void DrawDual::GVEntity::reset(const int index) 02220 { 02221 assert(index >= 0); 02222 dualSurfs[index] = 0; 02223 pointPos[index][0] = pointPos[index][1] = 0; 02224 vtkEntityIds[index] = -1; 02225 02226 int dim = MBI->dimension_from_handle(moabEntity); 02227 02228 // use gvizEdges to tell whether we're an edge or not 02229 if (0 == dim) { 02230 if (gvizPoints[index]) { 02231 if (useGraphviz) 02232 free(gvizPoints[index]); 02233 else 02234 MBI->delete_entities((EntityHandle*)&gvizPoints[index], 1); 02235 gvizPoints[index] = NULL; 02236 } 02237 } 02238 else if (1 == dim) { 02239 vtkEntityIds[index+2] = -1; 02240 vtkEntityIds[index+3] = -1; 02241 if (gvizEdges[index]) { 02242 if (useGraphviz) 02243 free(gvizEdges[index]); 02244 else { 02245 assert(gvizEdges[index]); 02246 MBI->delete_entities((EntityHandle*)&gvizEdges[index], 1); 02247 } 02248 gvizEdges[index] = NULL; 02249 } 02250 if (gvizEdges[index+2]) { 02251 if (useGraphviz) { 02252 free(gvizEdges[index+2]); 02253 if (gvizPoints[index]) free(gvizPoints[index]); 02254 } 02255 else { 02256 assert(gvizEdges[index+2]); 02257 MBI->delete_entities((EntityHandle*)&gvizEdges[index+2], 1); 02258 if (gvizPoints[index]) 02259 MBI->delete_entities((EntityHandle*)&gvizPoints[index], 1); 02260 } 02261 gvizEdges[index+2] = NULL; 02262 gvizPoints[index] = NULL; 02263 } 02264 } 02265 02266 myActors[index] = NULL; 02267 } 02268 02269 void DrawDual::get_graph_points(const EntityHandle *ents, const int num_ents, 02270 EntityHandle dual_surf, void **points) 02271 { 02272 std::vector<GVEntity *> gvs(num_ents); 02273 ErrorCode result = MBI->tag_get_data(gvEntityHandle, ents, num_ents, &gvs[0]); 02274 if (MB_SUCCESS != result) return; 02275 02276 for (int i = 0; i < num_ents; i++) { 02277 assert(NULL != gvs[i]); 02278 int index = gvs[i]->get_index(dual_surf); 02279 assert(index >= 0 && index < 3); 02280 points[i] = gvs[i]->gvizPoints[index]; 02281 } 02282 } 02283 02284 void DrawDual::get_graph_points(Range ents, 02285 EntityHandle dual_surf, void **points) 02286 { 02287 std::vector<GVEntity *> gvs(ents.size()); 02288 ErrorCode result = MBI->tag_get_data(gvEntityHandle, ents, &gvs[0]); 02289 if (MB_SUCCESS != result) return; 02290 02291 for (unsigned int i = 0; i < ents.size(); i++) { 02292 int index = gvs[i]->get_index(dual_surf); 02293 assert(index >= 0 && index < 3); 02294 points[i] = gvs[i]->gvizPoints[index]; 02295 } 02296 } 02297 02298 ErrorCode DrawDual::smooth_dual_surf(EntityHandle dual_surf, 02299 Range &dcells, 02300 Range &dedges, Range &dverts, 02301 Range &face_verts, Range &loop_edges) 02302 { 02303 // compute the starting positions of the boundary points, and fix them; 02304 // has to come after construct_graphviz_edges 'cuz we need edge gventities for 2-pt loops 02305 // copy face_verts 'cuz compute_fixed_points may add face verts for pillow sheets 02306 ErrorCode result = compute_fixed_points(dual_surf, dverts, 02307 face_verts, loop_edges); RR; 02308 02309 const int num_its = 10; 02310 02311 dverts = subtract( dverts, face_verts); 02312 double tmp_coords[12]; 02313 MeshTopoUtil mtu(vtkMOABUtils::mbImpl); 02314 02315 std::vector<EntityHandle> graph_points(dverts.size()); 02316 get_graph_points(dverts, dual_surf, (void**) &graph_points[0]); 02317 02318 // get any edge-adjacent points which are valence-2 and add to list of 02319 // smooth points; these are edge mid-pts, which should be smoothed so they 02320 // don't remain at the center 02321 for (unsigned int j = 0; j < dverts.size(); j++) { 02322 EntityHandle this_point = graph_points[j]; 02323 02324 // get all neighbor verts 02325 Range nverts, tmp_edges; 02326 result = mtu.get_bridge_adjacencies(this_point, 1, 0, nverts); RR; 02327 if (!(4 == nverts.size() || 3 == nverts.size())) { 02328 std::cerr << "Smoothing sheet failed; dumping file." << std::endl; 02329 vtkMOABUtils::dualTool->delete_whole_dual(); 02330 EntityHandle save_set; 02331 ErrorCode result = vtkMOABUtils::mbImpl->create_meshset(MESHSET_SET, save_set); 02332 if (MB_SUCCESS != result) return MB_FAILURE; 02333 Range hexes; 02334 result = vtkMOABUtils::mbImpl->get_entities_by_type(0, MBHEX, hexes); 02335 if (MB_SUCCESS != result) return MB_FAILURE; 02336 result = vtkMOABUtils::mbImpl->add_entities(save_set, hexes); 02337 if (MB_SUCCESS != result) return MB_FAILURE; 02338 vtkMOABUtils::mbImpl->write_file("tmp.h5m", NULL, NULL, &save_set, 1); 02339 return MB_FAILURE; 02340 } 02341 02342 for (Range::iterator rit = nverts.begin(); rit != nverts.end(); rit++) { 02343 tmp_edges.clear(); 02344 result = MBI->get_adjacencies(&(*rit), 1, 1, false, tmp_edges); RR; 02345 if (2 == tmp_edges.size() && 02346 std::find(graph_points.begin(), graph_points.end(), *rit) == graph_points.end()) 02347 graph_points.push_back(*rit); 02348 } 02349 } 02350 02351 std::vector<double> new_coords(3*graph_points.size()), 02352 old_coords(3*graph_points.size()); 02353 02354 for (int i = 0; i < num_its; i++) { 02355 // get starting coords for all verts 02356 if (0 == i) { 02357 result = MBI->get_coords(&graph_points[0], graph_points.size(), 02358 &old_coords[0]); RR; 02359 } 02360 else old_coords.swap(new_coords); 02361 02362 02363 for (unsigned int j = 0; j < graph_points.size(); j++) { 02364 EntityHandle this_point = graph_points[j]; 02365 02366 // get all neighbor verts 02367 Range nverts; 02368 result = mtu.get_bridge_adjacencies(this_point, 1, 0, nverts); RR; 02369 assert(4 == nverts.size() || 3 == nverts.size() || 2 == nverts.size()); 02370 02371 // get coords for those verts 02372 result = MBI->get_coords(nverts, &tmp_coords[0]); RR; 02373 02374 // compute new coords using inverse length-weighted laplacian 02375 double denom = 0.0, delta[3] = {0.0, 0.0, 0.0}; 02376 for (unsigned int k = 0; k < nverts.size(); k++) { 02377 double tdelta[3]; 02378 tdelta[0] = (tmp_coords[3*k] - old_coords[3*j]); 02379 tdelta[1] = tmp_coords[3*k+1] - old_coords[3*j+1]; 02380 tdelta[2] = tmp_coords[3*k+2] - old_coords[3*j+2]; 02381 double lsq = sqrt(tdelta[0]*tdelta[0] + tdelta[1]*tdelta[1] + 02382 tdelta[2]*tdelta[2]); 02383 if (true) lsq = 1.0; 02384 denom += lsq; 02385 for (int l = 0; l < 3; l++) delta[l] += lsq*tdelta[l]; 02386 } 02387 if (0 != denom) { 02388 for (int l = 0; l < 3; l++) 02389 new_coords[3*j+l] = old_coords[3*j+l] + delta[l]/denom; 02390 } 02391 else { 02392 for (int l = 0; l < 3; l++) 02393 new_coords[3*j+l] = old_coords[3*j+l]; 02394 } 02395 } 02396 02397 // set the new coordinate positions after each iteration 02398 result = MBI->set_coords(&graph_points[0], graph_points.size(), 02399 &new_coords[0]); RR; 02400 } 02401 02402 return MB_SUCCESS; 02403 } 02404 02405 ErrorCode DrawDual::process_pick(EntityHandle dual_surf, 02406 const double x, const double y, 02407 Range &picked_ents) 02408 { 02409 // get vertices on interior of 3d sheet 02410 Range int_verts, face_verts; 02411 ErrorCode result = vtkMOABUtils::dualTool->get_dual_entities(dual_surf, NULL, NULL, 02412 &int_verts, &face_verts, NULL); 02413 int_verts = subtract( int_verts, face_verts); 02414 02415 // get vertices on sheet drawing for that sheet, and their positions 02416 std::vector<EntityHandle> int_points(int_verts.size()); 02417 get_graph_points(int_verts, dual_surf, (void**) &int_points[0]); 02418 std::vector<double> point_coords(3*int_verts.size()); 02419 result = MBI->get_coords(&int_points[0], int_points.size(), 02420 &point_coords[0]); RR; 02421 02422 // get the closest point to those 02423 double dist, closest_coords[3]; 02424 EntityHandle closest_point; 02425 for (unsigned int i = 0; i < int_verts.size(); i++) { 02426 double this_dist = (point_coords[3*i] - x)*(point_coords[3*i] - x) + 02427 (point_coords[3*i+1] - y)*(point_coords[3*i+1] - y); 02428 if (0 == i || this_dist < dist) { 02429 dist = this_dist; 02430 closest_point = int_points[i]; 02431 for (int j = 0; j < 3; j++) closest_coords[j] = point_coords[3*i+j]; 02432 } 02433 } 02434 02435 // get connected edges, points, in star order 02436 std::vector<EntityHandle> conn_edges(4); 02437 EntityHandle dual_edges_3d[4], dual_curves[4], conn_verts[4]; 02438 result = MBI->get_adjacencies(&closest_point, 1, 1, false, conn_edges); RR; 02439 assert(4 == conn_edges.size()); 02440 // conn_edges are on sheet diag; each points to 3d dual edge through curve tag, 02441 // and each 3d dual edge points to dual curve set through that same tag 02442 result = MBI->tag_get_data(dualCurveTagHandle, &conn_edges[0], 4, 02443 dual_edges_3d); RR; 02444 result = MBI->tag_get_data(dualCurveTagHandle, dual_edges_3d, 4, 02445 dual_curves); RR; 02446 // should only be 2 distinct curves; don't handle all dual curves the same 02447 // for now, to avoid geometric check 02448 assert((dual_curves[0] == dual_curves[1] && dual_curves[0] != dual_curves[2]) || 02449 (dual_curves[0] == dual_curves[2] && dual_curves[0] != dual_curves[1]) || 02450 (dual_curves[0] == dual_curves[3] && dual_curves[0] != dual_curves[1])); 02451 02452 // if same curves are next to each other, switch edges 02453 if (dual_curves[0] == dual_curves[1]) { 02454 SWAP(conn_edges[1], conn_edges[2]); 02455 SWAP(dual_curves[1], dual_curves[2]); 02456 SWAP(dual_edges_3d[1], dual_edges_3d[2]); 02457 } 02458 02459 else if (dual_curves[1] == dual_curves[2]) { 02460 SWAP(conn_edges[0], conn_edges[1]); 02461 SWAP(dual_curves[0], dual_curves[1]); 02462 SWAP(dual_edges_3d[0], dual_edges_3d[1]); 02463 } 02464 assert(dual_curves[0] == dual_curves[2] && 02465 dual_curves[1] == dual_curves[3]); 02466 02467 // get connected points & their coords 02468 for (int i = 0; i < 4; i++) { 02469 const EntityHandle *connect; 02470 int num_connect; 02471 result = MBI->get_connectivity(conn_edges[i], connect, num_connect); RR; 02472 if (connect[0] == closest_point) conn_verts[i] = connect[1]; 02473 else if (connect[1] == closest_point) conn_verts[i] = connect[0]; 02474 else assert(false); 02475 } 02476 02477 double conn_coords[12]; 02478 result = MBI->get_coords(conn_verts, 4, conn_coords); RR; 02479 CartVect pick_closest(x - closest_coords[0], y - closest_coords[1], 0.0); 02480 CartVect conn_vects[4]; 02481 02482 // pre-compute normalized vectors from closest to each connected vertex 02483 for (int i = 0; i < 4; i++) { 02484 conn_vects[i] = CartVect(conn_coords[3*i]-closest_coords[0], 02485 conn_coords[3*i+1]-closest_coords[1], 0.0); 02486 conn_vects[i].normalize(); 02487 } 02488 02489 // now evaluate to see if we picked an edge or 2cell 02490 for (int i = 0; i < 4; i++) { 02491 double this_dot = conn_vects[i] % pick_closest; 02492 02493 // if positive dot prod & dist < some tol, picked the edge 02494 if (0.0 < this_dot && 02495 abs(dist - this_dot*this_dot) < 1) 02496 picked_ents.insert(dual_edges_3d[i]); 02497 02498 // else if cross products w/ neighboring vectors opposite, we're in the 2cell 02499 else if (0.0 < this_dot && 02500 (conn_vects[i] * pick_closest) % (conn_vects[(i+1)%4] * pick_closest) 02501 < 0.0) { 02502 Range common_2cell; 02503 EntityHandle edges[2] = {dual_edges_3d[i], dual_edges_3d[(i+1)%4]}; 02504 result = MBI->get_adjacencies(edges, 2, 2, false, common_2cell); 02505 assert(common_2cell.size() > 0); 02506 picked_ents.insert(*common_2cell.begin()); 02507 } 02508 } 02509 02510 return MB_SUCCESS; 02511 } 02512 02513 EntityHandle DrawDual::get_dual_surf(vtkRenderer *this_ren) 02514 { 02515 std::map<EntityHandle, GraphWindows>::iterator mit; 02516 02517 for (mit = surfDrawrings.begin(); mit != surfDrawrings.end(); mit++) { 02518 vtkRenderer *gw_ren = NULL; 02519 if ((*mit).second.sheetDiagram) 02520 gw_ren = (*mit).second.sheetDiagram->sheet_diagram()->GetRenderWindow()-> 02521 GetRenderers()->GetFirstRenderer(); 02522 if (gw_ren == this_ren) return (*mit).first; 02523 } 02524 02525 return 0; 02526 }