moab
DrawDual.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines