container.cc

Go to the documentation of this file.
00001 #include "container.h"
00002 
00003 #include <dirent.h>
00004 
00005 #include <cassert>
00006 #include <cmath>
00007 #include <fstream>
00008 #include <iostream>
00009 
00010 #include "box.h"
00011 #include "edge.h"
00012 #include "face.h"
00013 #include "gain_schedule.h"
00014 #include "intersecting_faces.h"
00015 #include "misc.h"
00016 #include "nice.h"
00017 #include "object.h"
00018 #include "point.h"
00019 #include "space.h"
00020 #include "state.h"
00021 #include "vertex.h"
00022 #include "virtual_disp.h"
00023 #include "log.h"
00024 
00025 using std::cout;
00026 using std::endl;
00027 
00028 Container * Container::only_one = NULL;
00029 
00030 Container & Container::instance()
00031 {
00032   // Not thread-safe.
00033   // -- lock mutex
00034   if (only_one == NULL)
00035     only_one = new Container();
00036   // -- unlock mutex
00037   return *only_one;
00038 }
00039 
00040 Container::Container (const Container& rhs)
00041   :fo(NULL),energy(0),force(0),files(),frozen(),o()
00042 {
00043   cout << "Copy constructor prohibited on instances of Container class.\n";
00044   cout << "Container " << rhs.energy << endl;
00045   exit(0);
00046 }
00047 
00048 Container& Container::operator= (const Container& rhs)
00049 {
00050   cout << "Copy assignment operator prohibited on instances of Container class.\n";
00051   cout << "Container " << rhs.energy << endl;
00052   exit(0);
00053 }
00054 
00055 Container::Container (void)
00056   :fo(NULL),energy(0),force(0),files(),frozen(),o()
00057 {
00058   scanDir();
00059   scanFiles();
00060   createEdges();
00061   findVertAdj();
00062   sortAdjacentFaces();
00063 }
00064 
00069 int Container::getVertexCount (void) const
00070 {
00071   int j = 0;
00072   // for each object
00073   for (o_cit i=o.begin();i!=o.end();i++)
00074   {
00075     j += i->v.size();
00076   }
00077   return j;
00078 }
00079 
00084 int Container::getFaceCount (void) const
00085 {
00086   int j = 0;
00087   // for each object
00088   for (o_cit i=o.begin();i!=o.end();i++)
00089   {
00090     j += i->f.size();
00091   }
00092   return j;
00093 }
00094 
00099 int Container::getEdgeCount (void) const
00100 {
00101   int j=0;
00102   // for each object
00103   for (o_cit i=o.begin();i!=o.end();i++)
00104   {
00105     j += i->e.size();
00106   }
00107   return j;
00108 }
00109 
00110 /* THIS WAS A DOXYGEN COMMENT
00111  * Write virtual displacement of each vertex in topN to file.
00112  *
00113  * \param[in] group Number of group of vertices being moved.
00114  * \param[in] suffix String to append to file name to identify
00115  * when in the program the file was written.
00116  *
00117  * \return Void.
00118  */
00119 
00120 //void Container::writeTopNSepDis (int group,str suffix)
00121 //{
00122 //  cout << "Write topN separation distances................";
00123 //  cout.flush();
00124 //  // open output file
00125 //  std::ofstream F;
00126 //  char filename[1024];
00127 //  sprintf(filename,"%stopN_sep_dis_%d_%s",OUTPUT_DATA_DIR.c_str(),group,suffix.c_str());
00128 //  F.open(filename);
00129 //  if(F.fail())
00130 //  {
00131 //    fprintf(stderr,"\nCouldn't open output file %s\n",filename);
00132 //    exit(0);
00133 //  }
00134 //  // for each topN
00135 //  for(tv_cit i=Virtual_Disp::instance().beginTopN();i!=Virtual_Disp::instance().endTopN();i++)
00136 //  {
00137 //    Vertex* vv = (*i).second;
00138 //    if(vv->getFace()!=NULL)
00139 //    {
00140 //      // get closest point
00141 //      //double pC[3];
00142 //      vec_d pC = getNearPtOnFaceToVertex(vv->getFace(),vv);
00143 //      // compute separation vector
00144 //      double ss[3];
00145 //      for(int j=0;j<3;j++){ ss[j]=pC[j]-vv->getCoord(j); }
00146 //      // compute separation distance //////////
00147 //      double sd = sqrt(dot(ss,ss));
00148 //      F << sd << endl;
00149 //    }
00150 //  }
00151 //  F.close();
00152 //  cout << "complete.\n";
00153 //  cout.flush();
00154 //}
00155 
00166 void Container::checkClosestFace (int group,str suffix)
00167 {
00168   cout << "Iteration " << group << ": Check closest face to vertices "
00169         << suffix << "...........";
00170   cout.flush();
00171   // open output file
00172   std::ofstream F;
00173   char filename[1024];
00174   sprintf(filename,"%sclosest_face_discrepencies_%d_%s",OUTPUT_DATA_DIR.c_str(),group,suffix.c_str());
00175   F.open(filename);
00176   if(F.fail())
00177   {
00178     fprintf(stderr,"\nCouldn't open output file %s\n",filename);
00179     exit(0);
00180   }
00181   // for each object in container
00182   for (o_it i=o.begin();i!=o.end();i++)
00183   {
00184     //if ((*i).getName()!="a108") continue;
00185     // for each vertex in object
00186     for (v_it j=i->v.begin();j!=i->v.end();j++)
00187     {
00188       //if ((*j).getIndex()!=1396) continue;
00189       if(binary_search(frozen.begin(),frozen.end(),&(*j))==false)
00190       {
00191         Face *ncl = NULL;
00192         vec_d p;
00193         double sqd=0.0;
00194         // find closest point to current vertex
00195         findClosestPtToVertex(&(*j),p,sqd,ncl);
00196         // if closest faces are the same, then no problem
00197         if(j->getFace()==ncl) continue;
00198         // closest face is incorrect
00199         cout << "Container::checkClosestFace: Error. "
00200               << "Stored closest face to vertex does not match actual closest face.\n";
00201         j->print(std::cout);
00202         if(ncl != NULL)
00203         {
00204           cout << "\nActual face (good)\n"
00205                 << "separation distance " << sqrt(sqd) << endl
00206                 << "closest point ["
00207                 << p[0] << " "
00208                 << p[1] << " "
00209                 << p[2] << "]\n";
00210           ncl->print(std::cout);
00211           cout << endl << endl;
00212         }
00213         else
00214         {
00215           cout << ", good NULL]\n\n";
00216         }
00217         if (j->getFace() != NULL)
00218         {
00219           vec_fp fp(1,const_cast<Face*>(j->getFace()));
00220           ncl=NULL;
00221           sqd=0.0;
00222           p[0]=p[1]=p[2]=0;
00223           if (findClosestPtToVertexAmongFaces(&(*j),fp,p,sqd,ncl)==true)
00224           {
00225             cout << "\nStored face (bad)\n"
00226                   << "separation distance " << sqrt(sqd) << endl
00227                   << "closest point ["
00228                   << p[0] << " "
00229                   << p[1] << " "
00230                   << p[2] << "]\n";
00231           }
00232           else
00233           {
00234             cout << "\nStored face (bad)\n"
00235                   << "no closest point was found on stored face." << endl;
00236           }
00237           j->getFace()->print(std::cout);
00238           cout << endl;
00239         }
00240         else
00241         {
00242           cout << "[bad NULL]\n";
00243         }
00244         writeMeshData(666);
00245         exit(0);
00246       }
00247     }
00248   }
00249   F.close();
00250   cout << "complete.\n";
00251   cout.flush();
00252 }
00253 
00260 void Container::createEdges (void)
00261 {
00262   if(!WRITE_VERBOSE_INIT)
00263   {
00264     cout << "Creating edges.................................";
00265     cout.flush();
00266   }
00267   int k=1;
00268   double goal = 0.2;
00269   printf("0%%..");
00270   fflush(stdout);
00271   // for each object, create edges
00272   double a = 1.0/o.size();
00273   for (o_it i=o.begin();i!=o.end();i++)
00274   {
00275     if(WRITE_VERBOSE_INIT)
00276     {
00277       cout << "Creating edges for " << i->getName() << "...";
00278       cout.flush();
00279     }
00280     i->e.reserve(3*i->v.size());
00281     // determine number of digits in largest vertex index
00282     int num_digits = i->setNumDigits();
00283     // create map for finding edges
00284     map_s_ep hm;
00285     // for each face
00286     for (f_it j=i->f.begin();j!=i->f.end();j++)
00287     {
00288       i->processEdge(&(*j),j->getVertex(0),j->getVertex(1),hm,num_digits);
00289       i->processEdge(&(*j),j->getVertex(1),j->getVertex(2),hm,num_digits);
00290       i->processEdge(&(*j),j->getVertex(2),j->getVertex(0),hm,num_digits);
00291     }
00292     if(WRITE_VERBOSE_INIT)
00293     {
00294       cout << "complete.\n";
00295       cout.flush();
00296     }
00297     // track progress
00298     double progress = static_cast<double>(k++)*a;
00299     if(progress>goal)
00300     {
00301       printf("%d%%..",static_cast<int>(goal*100));
00302       fflush(stdout);
00303       goal+=0.2;
00304     }
00305     assert(i->getFE()==&i->e[0]);
00306   }
00307   printf("100%%..");
00308   fflush(stdout);
00309   if(!WRITE_VERBOSE_INIT)
00310   {
00311     cout << "complete.\n";
00312     cout.flush();
00313   }
00314 }
00315 
00322 void Container::findVertAdj (void)
00323 {
00324   if(!WRITE_VERBOSE_INIT)
00325   {
00326     cout << "Finding vertex adjacencies.....................";
00327     cout.flush();
00328   }
00329   if (WRITE_NEIGHBORHOOD_TO_FILE)
00330   {
00331     // open dat file
00332     str log_file = OUTPUT_DATA_DIR + NEIGHBORHOOD_FILE;
00333     std::ofstream this_file (log_file.c_str());
00334     this_file.close();
00335   }
00336   // for each object, find vertex adjacencies
00337   if(WRITE_VERBOSE_INIT){ cout << endl;} 
00338   for (o_it i=o.begin();i!=o.end();i++)
00339   {
00340     if(WRITE_VERBOSE_INIT)
00341     {
00342       cout << "Finding adjacencies for " << i->getName() << "...";
00343       cout.flush();
00344     }
00345     i->findVertAdj();
00346     if(WRITE_VERBOSE_INIT)
00347     {
00348       cout << "complete.\n";
00349       cout.flush();
00350     }
00351   }
00352   if(!WRITE_VERBOSE_INIT)
00353   {
00354     cout << "complete.\n";
00355     cout.flush();
00356   }
00357 }
00358 
00364 void Container::scanDir (void)
00365 {
00366   str::size_type found;
00367   DIR *pdir;            // pointer to a directory data structure
00368   struct dirent *pent;  // pointer to dirent structure
00369   if(WRITE_VERBOSE_INIT){ cout << endl;cout.flush();}
00370   pdir = opendir(INPUT_DATA_DIR.c_str());
00371   if (!pdir) {printf("Error. Could not open %s.\n",INPUT_DATA_DIR.c_str());exit(1);}
00372   while ((pent=readdir(pdir)))
00373   {
00374     // copy char array to string
00375     str line = pent->d_name;
00376     // if file of typ *.mesh
00377     found = line.find(".mesh",0);
00378     // if found
00379     if (found != str::npos)
00380     {
00381       // save filename
00382       files.push_back(line);
00383       // update index
00384       // print file found to screen
00385       if(WRITE_VERBOSE_INIT)
00386       {
00387         cout << "file found: " << line << "\n"; cout.flush();
00388       }
00389     }
00390   }
00391   if (files.empty()==true)
00392   {
00393     cout << "\n\nContainer::scanDir: Error. No .mesh files were found in input dir.\n"
00394           << "Input directory = " << INPUT_DATA_DIR
00395           << endl << endl;
00396     exit(0);
00397   }
00398   closedir(pdir);
00399   if(WRITE_VERBOSE_INIT){ cout << endl;cout.flush();}
00400   o.reserve(files.size());
00401 }
00402 
00409 void Container::scanFiles (void)
00410 {
00411   cout << "\nInput Data Directory\n"
00412         << INPUT_DATA_DIR << "\n\n"
00413         << "Reading Directory..............................";
00414   cout.flush();
00415   // for each input file
00416   for (uint count=0;count<files.size();count++)
00417   {
00418     // record object name
00419     str::size_type pos1 = files[count].find(".",0);
00420     if (!(pos1==str::npos))
00421     {
00422         o.push_back(Object(files[count].substr(0,pos1)));
00423         if (count==0)
00424         {
00425           fo = &o[0];
00426         }
00427     }
00428     else
00429     {
00430       cout << "Error! Object name was not found in " << files[count] << "\n";
00431       exit(1);
00432     }
00433     // scan file
00434     str file = INPUT_DATA_DIR + files[count];
00435     if(WRITE_VERBOSE_INIT)
00436     {
00437       cout << "loading file " << files[count] << "...";
00438     }
00439     cout.flush();
00440     scanFile(&(o.back()),file.c_str());
00441     if(WRITE_VERBOSE_INIT)
00442     {
00443       cout << "complete.\n";
00444       cout.flush();
00445     }
00446   }
00447   // trim Object vector
00448   assert(fo==&o[0]);
00449   cout << "complete.\n";
00450 }
00451 
00458 void Container::assessFile (Object * const obj,const char *filename)
00459 {
00460   char line[2048];
00461   int v_count=0,f_count=0;
00462   // open file
00463   FILE *F = fopen(filename,"r");
00464   if (!F) { printf("Couldn't open input file %s\n",filename);return;}
00465   // for every line in file
00466   for (char *str=fgets(line,2048,F) ; str!=NULL ; str=fgets(line,2048,F))
00467   {
00468     // skip leading whitespace
00469     while (strchr(" \t,",*str)!=NULL) { str++;}
00470     // if first character is V for Vertex
00471     if (strchr("V",*str)!=NULL){ v_count++; }
00472     // if first character is F for Face
00473     else if (strchr("F",*str)!=NULL){ f_count++; }
00474   }
00475   fclose(F);
00476   obj->v.reserve(v_count);
00477   obj->f.reserve(f_count);
00478 }
00479 
00485 void Container::scanFile (Object * const obj,const char *filename)
00486 {
00487   char line[2048],*str;
00488   FILE *F;
00489   int vertex_num=0,polygon_num=0;
00490   // reserve space for vertices and faces
00491   assessFile(obj,filename);
00492   vec_vp vp;
00493   // open file
00494   F = fopen(filename,"r");
00495   if (!F) { printf("Couldn't open input file %s\n",filename);return;}
00496   // for every line in file
00497   for (str=fgets(line,2048,F) ; str!=NULL ; str=fgets(line,2048,F))
00498   {
00499     // skip leading whitespace
00500     while (strchr(" \t,",*str)!=NULL) { str++;}
00501     // if first character is V for Vertex,
00502     // add new linked list class instance
00503     if (strchr("V",*str)!=NULL)
00504     {
00505       obj->v.push_back(Vertex(str,obj));
00506       Vertex *vv = &obj->v.back();
00507       vp.push_back(vv);
00508       if(vertex_num==0)
00509       {
00510         obj->setFV(&obj->v[0]);
00511       }
00512       vertex_num++;
00513       if (PRINT_FLAG)
00514       {
00515         cout<<filename<<": ";
00516         vv->print(std::cout);
00517       }
00518     }
00519     // if first character is F for Face,
00520     // add new linked list class instance
00521     else if (strchr("F",*str)!=NULL)
00522     {
00523       obj->f.push_back(Face(str,vp));
00524       if(polygon_num==0)
00525       {
00526         obj->setFF(&obj->f[0]);
00527       }
00528       polygon_num++;
00529       if (PRINT_FLAG)
00530       {
00531         Face *ff = &obj->f.back();
00532         ff->print(std::cout);
00533       }
00534     }
00535   }
00536   assert(obj->getFV()==&obj->v[0]);
00537   assert(obj->getFF()==&obj->f[0]);
00538   fclose(F);
00539 }
00540 
00547 Object* Container::getObjectPointer (str name) const
00548 {
00549   // for each object* in container
00550   for(o_cit i=o.begin();i!=o.end();i++)
00551   {
00552     if(name==i->getName())
00553     {
00554       return const_cast<Object*>(&(*i));
00555     }
00556   }
00557   return NULL;
00558 }
00559 
00567 mmap_oi Container::loadFrozenMap(const char *filename,s_set & not_found)
00568 {
00569   mmap_oi frozen_map;
00570   // open file
00571   FILE *F = fopen(filename,"r");
00572   if (!F)
00573   {
00574     printf("Couldn't open input file %s\n",filename);
00575     return frozen_map;
00576   }
00577   // for every line in file
00578   // NOTE THE ASSUMED MAX character array length
00579   char val[2048],*eptr=NULL,line[2048],*input=NULL;
00580   Object *o1=NULL,*o2=NULL;
00581   for (input=fgets(line,2048,F) ; input!=NULL ; input=fgets(line,2048,F))
00582   {
00583     // get object name
00584     while (strchr(" \t,",*input)!=NULL) { input++; }
00585     int i=0;
00586     while (strchr(" \t",*input)==NULL)
00587     {
00588       val[i++] = *input++;
00589     }
00590     val[i]=0;
00591     // val contains object name
00592 
00593     // get Object*
00594     if(o2==NULL){o1 = getObjectPointer(val);}
00595     else
00596     {
00597       // if new object name same as old object name
00598       if(strcmp(o2->getName().c_str(),val)==false){o1=o2;}
00599       else {o1 = getObjectPointer(val);}
00600     }
00601     o2=o1;
00602     if (o1==NULL)
00603     {
00604       not_found.insert(val);
00605     }
00606     else
00607     {
00608       // get vertex index
00609       while (strchr(" \t,",*input)!=NULL) { input++; }
00610       i=0;
00611       while (strchr("0123456789",*input)!=NULL)
00612       {
00613         val[i++] = *input++;
00614       }
00615       val[i]=0;
00616       int vi = (int)strtod(val,&eptr);
00617       if (val==eptr)
00618       {
00619         vi=0;
00620         printf("Error in reading vertex index\n");
00621         return frozen_map;
00622       }
00623       // create map entry
00624       frozen_map.insert(std::make_pair(o1,vi));
00625     }
00626   }
00627   // close file
00628   fclose(F);
00629   return frozen_map;
00630 }
00631 
00636 void Container::readFrozen (const char *filename)
00637 {
00638   cout << "Read frozen vertices...........................";
00639   cout.flush();
00640   // load frozen map: object*->Vertex_index (int)
00641   oi_it front;
00642   // set of object names from frozen file not found as Object name
00643   s_set not_found;
00644   // read frozen vertex file 
00645   mmap_oi frozen_map = loadFrozenMap(filename,not_found);
00646 
00647   // for each element in multimap
00648   for(oi_it i=frozen_map.begin();i!=frozen_map.end();i++)
00649   {
00650     Object *oo=(*i).first;
00651     int t = (*i).second;
00652     int a = t-1;
00653     int b = oo->v[a].getIndex();
00654     do{
00655       if(t==b){frozen.push_back(&oo->v[a]);}
00656       else if (t<b){a++;}
00657       else if (t>b){a--;}
00658     } while(t!=b);
00659   }
00660   sort(frozen.begin(),frozen.end());
00661   cout << "complete.\n";
00662 
00663   // print objects not found
00664   if (not_found.empty()==false)
00665   {
00666     cout << "\nContainer::loadFrozenMap: Warning.\n"
00667          << "No matching Object* found in container "
00668          << "for following frozen objects:\n";
00669     for (ss_it j = not_found.begin();j!=not_found.end();j++)
00670     {
00671       cout << *j << endl;
00672     }
00673     cout << endl;
00674   }
00675 
00676   cout.flush();
00677 }
00678 
00683 void Container::writeMeshData (int group) const
00684 {
00685   // for each object
00686   for(o_cit i=o.begin();i!=o.end();i++)
00687   {
00688     // create output filename
00689     str file;
00690     if (APPEND_ITERATION_NUMBER_TO_MESH)
00691     {
00692       char sint[32];
00693       sprintf(sint,"_%i.mesh",group);
00694       str line(sint);
00695       file = OUTPUT_DATA_DIR + i->getName() + OUTPUT_SUFFIX + line;
00696     }
00697     else
00698     {
00699       str line(".mesh");
00700       file = OUTPUT_DATA_DIR + i->getName() + OUTPUT_SUFFIX + line;
00701     }
00702     // open output file
00703     std::ofstream newfile;
00704     newfile.open (file.c_str(),std::ios::trunc);
00705     if(newfile.is_open())
00706     {
00707       newfile.precision(12);
00708       // for each vertex in object
00709       for(v_cit j=i->v.begin();j!=i->v.end();j++)
00710       {
00711         // print index and final coordinates
00712         newfile << "Vertex "
00713               << j->getIndex() << " "
00714               << j->getCoord(0) << " "
00715               << j->getCoord(1) << " "
00716               << j->getCoord(2) << "\n";
00717       }
00718       // for each face in object
00719       for(f_cit k=i->f.begin();k!=i->f.end();k++)
00720       {
00721         newfile << "Face "
00722               << k->getIndex() << " "
00723               << k->getVertex(0)->getIndex() << " "
00724               << k->getVertex(1)->getIndex() << " "
00725               << k->getVertex(2)->getIndex() << endl;
00726       }
00727       newfile.close();
00728     }
00729     else
00730     {
00731       cout << "\nContainer::writeMeshData: Error:"
00732             << "Failed to open output file="
00733             << file << endl;
00734       cout.flush();
00735     }
00736   }
00737 }
00738 
00744 void Container::checkFacesInBoxes (int group, str suffix)
00745 {
00746   cout << "Iteration " << group << ": Check assignment of faces to boxes "
00747         << suffix << ".......";
00748   cout.flush();
00749   vec_bp bp;
00750   // for each object in container
00751   for (o_it i=o.begin();i!=o.end();i++)
00752   {
00753     // for each face in object
00754     for (f_it j=i->f.begin();j!=i->f.end();j++)
00755     {
00756       // identify boxes that overlap to any degree the face bounding box.
00757       // This conservative approach, i.e. always including the actual 
00758       // face-intersecting boxes plus others is meant to save time,
00759       // Assumes the time to exclude the other boxes is greater than the
00760       // time to check intersection with other boxes later.
00761       bp.clear();
00762       Space::instance().getBoxesToCheck(&(*j),bp);
00763       // check  
00764       if (bp.empty())
00765       { 
00766         cout << "Container::checkFacesInBoxes: ERROR: NO BOXES FOR\n";
00767         j->print(std::cout);
00768         exit(1);
00769       }
00770       // check each box in bp to ensure Face* is present
00771       for(bp_it k=bp.begin();k!=bp.end();k++)
00772       {
00773         bool found = false;
00774         // find this face in this box
00775         for(fp_cit m=(*k)->begin();m!=(*k)->end();m++)
00776         {
00777           if((*m)==&(*j)){found=true;break;}
00778         }
00779         // if not found, then error
00780         if(found==false)
00781         {
00782           cout << "\n\nContainer::checkFacesInBoxes: "
00783                 << "Error. Face not found in box!\n";
00784           exit(0);
00785 
00786         }
00787       }
00788       // Note this will not detect erroneous Face*s in Boxes
00789       // but this is not as bad a problem as missing Face*s in boxes
00790     }
00791   }
00792   cout << "complete.\n";
00793   cout.flush();
00794 }
00795 
00800 void Container::assignFacesToBoxes (void)
00801 {
00802   cout << "Assign faces to boxes..........................";
00803   cout.flush();
00804   vec_bp bp;
00805   // clear boxes
00806   Space::instance().clearBoxes();
00808   // for each object in container
00809   for (o_it i=o.begin();i!=o.end();i++)
00810   {
00811     // for each face in object
00812     for (f_it j=i->f.begin();j!=i->f.end();j++)
00813     {
00814       // identify boxes that overlap to any degree the face bounding box.
00815       // This conservative approach, i.e. always including the actual 
00816       // face-intersecting boxes plus others is meant to save time,
00817       // assuming the time to exclude the other boxes is greater than the
00818       // time to check intersection with other boxes later.
00819       bp.clear();
00820       Space::instance().getBoxesToCheck(&(*j),bp);
00821       // check  
00822       if (bp.empty())
00823       { 
00824         cout << "ERROR: NO BOXES FOR\n";
00825         j->print(std::cout);
00826         exit(1);
00827       }
00828       // record face in boxes class
00829       Space::instance().recordFace(bp,&(*j));
00830       // update statistics on number of boxes per face
00831       Log::instance().updatePartitioningFaceStats(bp.size());
00832     }
00833   }
00834   cout << "complete.\n";
00835   cout.flush();
00836 }
00837 
00841 void Container::findClosestFaceToEachVertex (void)
00842 {
00843   cout << "Find closest face to each vertex...............";
00844   cout.flush();
00845   int index=1;
00846   double goal = 0.2;
00847   printf("0%%..");
00848   fflush(stdout);
00849   // for each object in container
00850   double a = 1.0/o.size();
00851   for (o_it i=o.begin();i!=o.end();i++)
00852   {
00853     // for each vertex in object
00854     for (v_it j=i->v.begin();j!=i->v.end();j++)
00855     {
00856       if (vertexIsFrozen(&(*j))==false)
00857       {
00859         // true => record closest face in vertex class
00860         // dummy => return closest face
00861         vec_d dummy_p;
00862         double dummy_sqd=0.0;
00863         Face *cl=NULL;
00864         if(findClosestPtToVertex(&(*j),dummy_p,dummy_sqd,cl)==true)
00865         {
00866           j->setFace(cl);
00867         }
00868         assert(dummy_sqd<1E10);
00869         assert(dummy_sqd==dummy_sqd);
00870       }
00871     }
00872     // track progress
00873     double progress = static_cast<double>(index++)*a;
00874     if(progress>goal)
00875     {
00876       printf("%d%%..",static_cast<int>(goal*100));
00877       fflush(stdout);
00878       goal+=0.2;
00879     }
00880   }
00881   printf("100%%..");
00882   fflush(stdout);
00883   cout << "complete.\n";
00884   cout.flush();
00885 }
00886 
00894 void Container::getFacesToCheckForClosestPt (vec_bp &bp,
00895                                              Vertex const * const v,
00896                                              vec_fp &cf) const
00897 {
00898   cf.clear();
00899   // for each box in search
00900   for (bp_it j=bp.begin();j!=bp.end();j++)
00901   {
00902     // for each face in box
00903     for (fp_cit k=(*j)->begin();k!=(*j)->end();k++)
00904     {
00905       // if face is not adjacent to current vertex
00906       // and not found in vector of old face pointers, of
00907       //if (v->faceIsAdjacent(*k)==false
00908       //&& binary_search(of.begin(),of.end(),*k)==false)
00909       if (v->faceIsAdjacent(*k)==false)
00910       {
00911         cf.push_back(*k);
00912       }
00913     }
00914   }
00915 }
00916 
00925 void Container::getSqDistFromVertexToSearchLimits (double * const & dist_squared,
00926                                                    double const * const search_region,
00927                                                    Vertex const * const v) const
00928 {
00929   // [xmin xmax ymin ymax zmin zmax]
00930   for (int i=0;i<6;i++)
00931   {
00932     int j = static_cast<int>(floor(static_cast<float>(i)/2.0));
00933     double a = search_region[i]-v->getCoord(j);
00934     dist_squared[i] = a*a;
00935   }
00936 }
00937 
00944 vec_d Container::getSearchRegion (Vertex const * const v) const
00945 {
00946   // vertex normal
00947   vec_d n = v->getNormal();
00948   // set normal length equal to SEARCH_RADIUS 
00949   double l = sqrt(dot(n,n));
00950   double m = sqrt(SEARCH_RADIUS_SQ)/l;
00951   n[0] = n[0]*m;
00952   n[1] = n[1]*m;
00953   n[2] = n[2]*m;
00954   // DEBUG
00955   //cout << "\nContainer::getSearchRegion: normal vector = ["
00956   //      << n[0] << " " 
00957   //      << n[1] << " " 
00958   //      << n[2] << "]\n"; 
00959   //cout << "\nContainer::getSearchRegion: normal vector length = "
00960   //<< l << endl;
00961   //cout << "\nContainer::getSearchRegion: normal vector scaling = "
00962   //<< m << endl;
00963   // DEBUG
00964   // length of rotation vectors
00965   double a = n[0]*n[0];
00966   double b = n[1]*n[1];
00967   double c = n[2]*n[2];
00968   double lrx = sqrt(b+c);
00969   double lry = sqrt(a+c);
00970   double lrz = sqrt(a+b);
00971   // DEBUG
00972   //cout << "\nContainer::getSearchRegion: rotation vector lengths = ["
00973   //      << lrx << " "
00974   //      << lry << " "
00975   //      << lrz << "]\n";
00976   // DEBUG
00977   // reducing number of product operations
00978   double sym[3] = {v->getCoord(0)+CLOSEST_POINT_COSINE*n[0], 
00979                    v->getCoord(1)+CLOSEST_POINT_COSINE*n[1], 
00980                    v->getCoord(2)+CLOSEST_POINT_COSINE*n[2]}; 
00981   double bol[3] = {n[0]*n[1],n[0]*n[2],n[1]*n[2]};
00982   // points for counterclockwise rotation
00983   vec_d myccwX(3,0.0);
00984   double d = CLOSEST_POINT_SINE/lrx;
00985   myccwX[0] = sym[0]+d*(lrx*lrx); 
00986   myccwX[1] = sym[1]+d*(-bol[0]); 
00987   myccwX[2] = sym[2]+d*(-bol[1]); 
00988   vec_d myccwY(3,0.0);
00989   double e = CLOSEST_POINT_SINE/lry;
00990   myccwY[0] = sym[0]+e*(-bol[0]); 
00991   myccwY[1] = sym[1]+e*(lry*lry); 
00992   myccwY[2] = sym[2]+e*(-bol[2]); 
00993   vec_d myccwZ(3,0.0);
00994   double f = CLOSEST_POINT_SINE/lrz;
00995   myccwZ[0] = sym[0]+f*(-bol[1]); 
00996   myccwZ[1] = sym[1]+f*(-bol[2]); 
00997   myccwZ[2] = sym[2]+f*(lrz*lrz); 
00998   // DEBUG
00999   //cout << "\nContainer::getSearchRegion: counterclockwise pt = ["
01000   //      << myccwX[0] << " "
01001   //      << myccwX[1] << " "
01002   //      << myccwX[2] << "]\n";
01003   //cout << "\nContainer::getSearchRegion: counterclockwise pt = ["
01004   //      << myccwY[0] << " "
01005   //      << myccwY[1] << " "
01006   //      << myccwY[2] << "]\n";
01007   //cout << "\nContainer::getSearchRegion: counterclockwise pt = ["
01008   //      << myccwZ[0] << " "
01009   //      << myccwZ[1] << " "
01010   //      << myccwZ[2] << "]\n";
01011   // DEBUG
01012   // points for clockwise rotation
01013   vec_d mycwX(3,0.0);
01014   double g = -CLOSEST_POINT_SINE/lrx;
01015   mycwX[0] = sym[0]+g*(lrx*lrx); 
01016   mycwX[1] = sym[1]+g*(-bol[0]); 
01017   mycwX[2] = sym[2]+g*(-bol[1]); 
01018   vec_d mycwY(3,0.0);
01019   double h = -CLOSEST_POINT_SINE/lry;
01020   mycwY[0] = sym[0]+h*(-bol[0]); 
01021   mycwY[1] = sym[1]+h*(lry*lry); 
01022   mycwY[2] = sym[2]+h*(-bol[2]); 
01023   vec_d mycwZ(3,0.0);
01024   double i = -CLOSEST_POINT_SINE/lrz;
01025   mycwZ[0] = sym[0]+i*(-bol[1]); 
01026   mycwZ[1] = sym[1]+i*(-bol[2]); 
01027   mycwZ[2] = sym[2]+i*(lrz*lrz); 
01028   // DEBUG
01029   //cout << "\nContainer::getSearchRegion: clockwise pt = ["
01030   //      << mycwX[0] << " "
01031   //      << mycwX[1] << " "
01032   //      << mycwX[2] << "]\n";
01033   //cout << "\nContainer::getSearchRegion: clockwise pt = ["
01034   //      << mycwY[0] << " "
01035   //      << mycwY[1] << " "
01036   //      << mycwY[2] << "]\n";
01037   //cout << "\nContainer::getSearchRegion: clockwise pt = ["
01038   //      << mycwZ[0] << " "
01039   //      << mycwZ[1] << " "
01040   //      << mycwZ[2] << "]\n";
01041   // DEBUG
01042   // find min and max values for x, y, z
01043   vec_d xval(7,0.0);
01044   xval[0] =      myccwX[0];
01045   xval[1] =      myccwY[0];
01046   xval[2] =      myccwZ[0];
01047   xval[3] =       mycwX[0];
01048   xval[4] =       mycwY[0];
01049   xval[5] =       mycwZ[0];
01050   xval[6] = v->getCoord(0);
01051   vec_d yval(7,0.0);
01052   yval[0] =      myccwX[1];
01053   yval[1] =      myccwY[1];
01054   yval[2] =      myccwZ[1];
01055   yval[3] =       mycwX[1];
01056   yval[4] =       mycwY[1];
01057   yval[5] =       mycwZ[1];
01058   yval[6] = v->getCoord(1);
01059   vec_d zval(7,0.0);
01060   zval[0] =      myccwX[2];
01061   zval[1] =      myccwY[2];
01062   zval[2] =      myccwZ[2];
01063   zval[3] =       mycwX[2];
01064   zval[4] =       mycwY[2];
01065   zval[5] =       mycwZ[2];
01066   zval[6] = v->getCoord(2);
01067   sort(xval.begin(),xval.end());
01068   sort(yval.begin(),yval.end());
01069   sort(zval.begin(),zval.end());
01070   // build search region
01071   vec_d sr(6,0.0);
01072   sr[0] = xval[0];
01073   sr[1] = xval[6];
01074   sr[2] = yval[0];
01075   sr[3] = yval[6];
01076   sr[4] = zval[0];
01077   sr[5] = zval[6];
01078   return sr;
01079 }
01080 
01090 bool Container::findClosestPtToVertex (Vertex const * const v,
01091                                        vec_d & p,
01092                                        double & sqd,
01093                                        Face *&ncl) const
01094 {
01095   Space & s(Space::instance());
01096   // get limits of search region using search cone
01097   // based on current vertex and vertex normal
01098   vec_d sr = getSearchRegion(v);
01099   // DEBUG
01100   //cout << "\nContainer::findClosestPtToVertex: search region = ["
01101   //      << sr[0] << " "
01102   //      << sr[1] << " "
01103   //      << sr[2] << " "
01104   //      << sr[3] << " "
01105   //      << sr[4] << " "
01106   //      << sr[5] << "]\n";
01107   // DEBUG
01108   // get box indices for search regeion : [xmin xmax ymin ymax zmin zmax]
01109   int indices[6] = { s.screenIndex(s.location2Index(sr[0],"x"),"x"),
01110                      s.screenIndex(s.location2Index(sr[1],"x"),"x"),
01111                      s.screenIndex(s.location2Index(sr[2],"y"),"y"),
01112                      s.screenIndex(s.location2Index(sr[3],"y"),"y"),
01113                      s.screenIndex(s.location2Index(sr[4],"z"),"z"),
01114                      s.screenIndex(s.location2Index(sr[5],"z"),"z")};
01115   // DEBUG
01116   //cout << "Container::findClosestPtToVertex: box indices = ["
01117   //      << indices[0] << " "
01118   //      << indices[1] << " "
01119   //      << indices[2] << " "
01120   //      << indices[3] << " "
01121   //      << indices[4] << " "
01122   //      << indices[5] << "]\n";
01123   // DEBUG
01124   // collect Box*s from indices
01125   vec_bp nb;
01126   nb.reserve(VECTOR_RESERVE);
01127   Space::instance().getBoxesFromIndices(NULL,indices,nb);
01128   // DEBUG
01129   //cout << "Container::findClosestPtToVertex: number of boxes = " << nb.size() << endl;
01130   // DEBUG
01131   // collect face*s vectors from boxes
01132   vec_fp cf;
01133   cf.reserve(VECTOR_RESERVE);
01134   getFacesToCheckForClosestPt(nb,v,cf);
01135   // sort and keep unique Face*
01136   sort(cf.begin(),cf.end());
01137   cf.assign(cf.begin(),unique(cf.begin(),cf.end()));
01138   // DEBUG
01139   //cout << "Container::findClosestPtToVertex: number of faces = " << cf.size() << endl;
01140   // DEBUG
01141   // count of face and boxes used in search
01142   Log::instance().updateClosestPtStats(cf.size(),nb.size());
01143   // search for closest point (cp) and closest distance (cd) among faces
01144   sqd=0.0;
01145   return findClosestPtToVertexAmongFaces(v,cf,p,sqd,ncl);
01146 }
01147 
01148 /* THIS WAS A DOXYGEN COMMENT
01149  * Find closest point to vertex among input faces.
01150  *
01151  * \param[in] v Vertex of interest.
01152  * \param[in] faces Faces on which closest point to vertex will be identified.
01153  * \param[out] p Closest point position, if found.
01154  * \param[out] sqd Squared distance between vertex
01155  * and closest point, if found.
01156  * \param[out] ncl Parent face of closest point, if found.
01157  * \return True if closest point found; otherwise false.
01158  */
01159 
01160 bool Container::findClosestPtToVertexAmongFaces (Vertex const * const v,
01161                                                  vec_fp & faces,
01162                                                  vec_d & p,
01163                                                  double & squareD,
01164                                                  Face * & ncl) const
01165 {
01166   bool gate = false;
01167   // if candidate faces were NOT found
01168   if (faces.empty()==true) return gate;
01169   // get vertex normal
01170   vec_d vn = v->getNormal();
01171   // collect vertices, edges
01172   vec_vp verts;
01173   vec_ep edges;
01174   // for each candidate face
01175   for (fp_it j=faces.begin();j!=faces.end();j++)
01176   {
01177     assert(*j!=NULL);
01178     // collect vertices and edges from candidate faces
01179     verts.push_back((*j)->getVertex(0));
01180     verts.push_back((*j)->getVertex(1));
01181     verts.push_back((*j)->getVertex(2));
01182     edges.push_back((*j)->getEdge(0));
01183     edges.push_back((*j)->getEdge(1));
01184     edges.push_back((*j)->getEdge(2));
01185   }
01186   map_df fmap;
01187   map_dt tmap;
01188   getDistFromVertexToVertices(v,verts,fmap,tmap);
01189   getDistFromVertexToEdges(v,edges,fmap,tmap);
01190   getDistFromVertexToFaces(v,faces,fmap,tmap);
01191   // invert vertex normal if vertex is not nice
01192   double vn_copy[3];
01193   if (Nice::instance().vertexIsNice(v))
01194   {
01195     for (int i=0;i<3;i++) vn_copy[i]=vn[i];
01196   } 
01197   else
01198   {
01199     for (int i=0;i<3;i++) vn_copy[i]=-vn[i];
01200   } 
01201   double a = dot(vn_copy,vn_copy);
01202   // MAKE SURE MAP IS SORTED SMALLEST TO LARGEST
01203   // for each triplet in map
01204   mdt_it i = tmap.begin();
01205   mdf_it j = fmap.begin();
01206   mdf_it k = j;
01207   while(i!=tmap.end())
01208   {
01209     // closest point
01210     //double a = (*i).second[0];
01211     //double b = (*i).second[1];
01212     //double c = (*i).second[2];
01213     p = (*i).second;
01214     // compute separation vector
01215     double sep_vec[3];
01216     sep_vec[0] = p[0]-v->getCoord(0);
01217     sep_vec[1] = p[1]-v->getCoord(1);
01218     sep_vec[2] = p[2]-v->getCoord(2);
01219     // compute cosine of angle between outward normal and separation vector
01220     // which is equal to dot product of vectors divided by vector magnitudes
01221     double sign = dot(sep_vec,vn_copy);
01222     if (sign>0.0)
01223     {
01224       double cos_angle2 = sign*sign/(a*dot(sep_vec,sep_vec));
01225       // is closest point located within search cone
01226       // with angle window as defined in controls.cc?
01227       // OR located very close to current vertex?
01228       if ( (cos_angle2 > CLOSEST_POINT_COSINE*CLOSEST_POINT_COSINE)
01229            || ((*i).first<(1.0*SCALE)))
01230       {
01231         // if square of separation distance is less than max allowed
01232         if ((*i).first<(SEARCH_RADIUS_SQ*SCALE*SCALE))
01233         {
01234           // if square of separation distance is less than
01235           // previously saved square of separation distance
01236           if ((*i).first<squareD||!squareD)
01237           {
01238             // DEBUG
01239             //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true)
01240             //{
01241             //  cout << "\nA closest point was found a distance ("
01242             //        << sqrt((*i).first)
01243             //        << ") away from target vertex on following face.\n";
01244             //  (*j).second->print(cout);
01245             //}
01246             // DEBUG
01247             // save 
01248             //if (update==true){ v->setFace((*j).second);}
01249             ncl=(*j).second;
01250             squareD = (*i).first;
01251             gate=true;
01252             break;
01253           }
01254         }
01255       }
01256     }
01257     i++;
01258     j++;
01259   }
01260   return gate;
01261 }
01262 
01272 bool Container::checkNearPtOnFaceInteriorToVertex (Face const * const f,
01273                                                   Vertex const * const v,
01274                                                   vec_d const & n) const
01275 {
01276   // describe line through vertex with same slope as n
01277   vec_d origin(3,0.0);
01278   origin[0] = v->getCoord(0);
01279   origin[1] = v->getCoord(1);
01280   origin[2] = v->getCoord(2);
01281   vec_d end(3,0.0);
01282   end[0]    = v->getCoord(0)+n[0];
01283   end[1]    = v->getCoord(1)+n[1];
01284   end[2]    = v->getCoord(2)+n[2];
01285   // determine if line intersects face f
01286   // last argument == true, ignore line intersection flag
01287   // ie, treat origin and end as line, not a line segment
01288   result r = checkLineFaceInt(f,origin,end,true);
01289   return r.poly_flag;
01290 }
01291 
01299 void Container::getNearPtOnFaceEdgesToVertex (Vertex const * const v,
01300                                               Face const * const f,
01301                                               Point & p) const
01302 {
01303   vec_d q;
01304   for (int i=0;i<3;i++)
01305   {
01306     if(const_cast<Face*>(f)->getEdge(i)->getNearPtToVertex(v,q)==true)
01307     {
01308       p.add(q[0],q[1],q[2]);
01309     }
01310   }
01311 }
01312 
01313 void Container::getDistFromVertexToVertices (Vertex const * const v,vec_vp & verts,map_df & fmap,map_dt & tmap) const
01314 {
01315   // sort and keep unique vertices
01316   sort(verts.begin(),verts.end());
01317   verts.assign(verts.begin(),unique(verts.begin(),verts.end()));
01318   // for each vertex
01319   for(vp_it i=verts.begin();i!=verts.end();i++)
01320   {
01321     // squared distance between v and *i
01322     //double sqd = ((*i)->getCoord(0)-v->getCoord(0))*((*i)->getCoord(0)-v->getCoord(0))+
01323     //             ((*i)->getCoord(1)-v->getCoord(1))*((*i)->getCoord(1)-v->getCoord(1))+
01324     //             ((*i)->getCoord(2)-v->getCoord(2))*((*i)->getCoord(2)-v->getCoord(2));
01325     double a = (*i)->getCoord(0)-v->getCoord(0);
01326     double b = (*i)->getCoord(1)-v->getCoord(1);
01327     double c = (*i)->getCoord(2)-v->getCoord(2);
01328     double sqd = a*a+b*b+c*c;
01329     // add squared distance to fmap
01330     fmap[sqd]=(*i)->getFaceNotAdjToVertex(v);
01331     // add vertex to tmap
01332     //vec_d t;
01333     //t.push_back((*i)->getCoord(0));
01334     //t.push_back((*i)->getCoord(1));
01335     //t.push_back((*i)->getCoord(2));
01336     //if(!tmap.count(sqd)){tmap[sqd]=t;}
01337     tmap[sqd]=*((*i)->getPos());
01338     // DEBUG
01339     //if(   v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01340     //   ((*i)->isMatch(TARGET_VERTEX_INDEX_2,TARGET_VERTEX_NAME_2)==true)
01341     //   )
01342     //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true)
01343     //{
01344     //  cout << "\nVertex " << (*i)->getIndex()
01345     //        << " searched for closest point to target vertex"
01346     //        << " and located a distance ("
01347     //        << sqrt(sqd) << ") away and associated with the following face.\n";
01348     //  (*i)->getFaceNotAdjToVertex(v)->print(cout);
01349     //}
01350     // DEBUG
01351   }
01352 }
01353 
01354 void Container::getDistFromVertexToEdges (Vertex const * const v,vec_ep & edges,map_df & fmap,map_dt & tmap) const
01355 {
01356   // sort and keep unique edges
01357   sort(edges.begin(),edges.end());
01358   edges.assign(edges.begin(),unique(edges.begin(),edges.end()));
01359   // for each edge
01360   for(ep_it i=edges.begin();i!=edges.end();i++)
01361   {
01362     // DEBUG
01363     //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01364     //   (
01365     //    (
01366     //    (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01367     //    (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01368     //    ) ||
01369     //    (
01370     //    (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01371     //    (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01372     //    )
01373     //   )
01374     //   )
01375     //{
01376     //  cout << "\nTarget edge searched for closest point to target vertex.\n";
01377     //}
01378     // DEBUG
01379     vec_d p;
01380     if((*i)->getNearPtToVertex(v,p)==true)
01381     {
01382       // squared distance between v and i
01383       double a = p[0]-v->getCoord(0);
01384       double b = p[1]-v->getCoord(1);
01385       double c = p[2]-v->getCoord(2);
01386       double sqd = a*a+b*b+c*c;
01387       // add squared distance to fmap
01388       fmap[sqd]=(*i)->getFaceNotAdjToVertex(v);
01389       // add vertex to tmap
01390       //vec_d t(3,0.0);
01391       //t[0] = ix;
01392       //t[1] = iy;
01393       //t[2] = iz;
01394       //if(!tmap.count(sqd)){tmap[sqd]=t;}
01395       tmap[sqd]=p;
01396       // DEBUG
01397       //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01398       //   (
01399       //    (
01400       //    (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01401       //    (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01402       //    ) ||
01403       //    (
01404       //    (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01405       //    (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01406       //    )
01407       //   )
01408       //   )
01409       //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true)
01410       //{
01411       //  cout << "Point ["
01412       //        << p[0] << " "
01413       //        << p[1] << " "
01414       //        << p[2] << "] found on edge "
01415       //        << " a distance ("
01416       //        << sqrt(sqd) << ") away and associated to following face\n";
01417       //  (*i)->getFaceNotAdjToVertex(v)->print(cout);
01418       //}
01419       // DEBUG
01420     }
01421     //else
01422     //{
01423     //  // DEBUG
01424     //  if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01425     //     (
01426     //      (
01427     //      (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01428     //      (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01429     //      ) ||
01430     //      (
01431     //      (*i)->getFace2()->getIndex()==TARGET_EDGE_FACE_INDEX_1 &&
01432     //      (*i)->getFace1()->getIndex()==TARGET_EDGE_FACE_INDEX_2 
01433     //      )
01434     //     )
01435     //     )
01436     //  {
01437     //    cout << "Point NOT found on edge.\n";
01438     //  }
01439     //  // DEBUG
01440     //}
01441   }
01442 }
01443 
01444 void Container::getDistFromVertexToFaces (Vertex const * const v,vec_fp & faces,map_df & fmap,map_dt & tmap) const
01445 {
01446   // sort and keep unique faces
01447   sort(faces.begin(),faces.end());
01448   faces.assign(faces.begin(),unique(faces.begin(),faces.end()));
01449   // initialize point class instance with current vertex
01450   Point p(v->getCoord(0),v->getCoord(1),v->getCoord(2));
01451   // for each face
01452   for(fp_it i=faces.begin();i!=faces.end();i++)
01453   {
01454     // DEBUG
01455     //if(   v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01456     //      (
01457     //   ((*i)->isMatch(TARGET_FACE_INDEX_1,TARGET_FACE_NAME)==true) ||
01458     //   ((*i)->isMatch(TARGET_FACE_INDEX_2,TARGET_FACE_NAME)==true)
01459     //   )
01460     //   )
01461     //{
01462     //  cout << "\nTarget face " << (*i)->getIndex()
01463     //        << " searched in full search for closest point to target vertex.\n";
01464     //}
01465     // DEBUG
01466     // get face normal
01467     vec_d fn = (*i)->getNormal();
01468     // compute vector connecting arbitrary face vertex and current vertex
01469     vec_d diff(3,0);
01470     diff[0] = (*i)->getVertex(0)->getCoord(0)-v->getCoord(0);
01471     diff[1] = (*i)->getVertex(0)->getCoord(1)-v->getCoord(1);
01472     diff[2] = (*i)->getVertex(0)->getCoord(2)-v->getCoord(2);
01473     // compute indicators
01474     double num=dot(fn,diff);
01475     // if current vertex does not lie on face plane
01476     bool inside = false;
01477     p.clear();
01478     if (fabs(num)>DOUBLE_EPSILON)
01479     {
01480       //compute point on face plane that is closest to current vertex
01481       // i.e. intersection of plane with face normal through current vertex
01482       // add to p if intersection is on face
01483       // return true if point lies inside or on an edge of face, false otherwise
01484       inside = checkNearPtOnFaceInteriorToVertex(*i,v,fn);
01485     }
01486     // if a point has been found, i.e. inside==true
01487     if(inside==true)
01488     {
01489       double uu=num/dot(fn,fn);
01490       p.add(v->getCoord(0)+uu*fn[0],v->getCoord(1)+uu*fn[1],v->getCoord(2)+uu*fn[2]);
01491       vec_d t;
01492       t.push_back(p.getClosestX());
01493       t.push_back(p.getClosestY());
01494       t.push_back(p.getClosestZ());
01495       // DEBUG
01496       //if(   v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01497       //      (
01498       //   ((*i)->isMatch(TARGET_FACE_INDEX_1,TARGET_FACE_NAME)==true) ||
01499       //   ((*i)->isMatch(TARGET_FACE_INDEX_2,TARGET_FACE_NAME)==true)
01500       //   )
01501       //   )
01502       //if(v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true)
01503       //{
01504       //  cout << "Point ("
01505       //        << t[0] << " "
01506       //        << t[1] << " "
01507       //        << t[2] << ") found on target face a distance ("
01508       //        << sqrt(p.getSqD()) << ") from target vertex.\n";
01509       //}
01510       // DEBUG
01511       // add squared distance to fmap
01512       fmap[p.getSqD()]=*i;
01513       // add point to tmap
01514       if(!tmap.count(p.getSqD())){tmap[p.getSqD()]=t;}
01515     }
01516     //else
01517     //{
01518     //  // DEBUG
01519     //  if(   v->isMatch(TARGET_VERTEX_INDEX_1,TARGET_VERTEX_NAME_1)==true &&
01520     //        (
01521     //     ((*i)->isMatch(TARGET_FACE_INDEX_1,TARGET_FACE_NAME)==true) ||
01522     //     ((*i)->isMatch(TARGET_FACE_INDEX_2,TARGET_FACE_NAME)==true)
01523     //     )
01524     //     )
01525     //  {
01526     //    cout << "Point not found on target face.\n";
01527     //  }
01528     //  // DEBUG
01529     //}
01530   }
01531 }
01532 
01533 double Container::getMinEdgeAngle (void) const
01534 {
01535   //cout << "Check edge angles..............................";
01536   //cout.flush();
01537   double min = 1E30;
01538   // for each object* in container
01539   for (o_cit i=o.begin();i!=o.end();i++)
01540   {
01541     // for each Edge* in object
01542     for (e_cit j=i->e.begin();j!=i->e.end();j++)
01543     {
01544       // record if angle is smallest so far
01545       double a = j->getAngle(); 
01546       if (a<min) min=a;
01547     }
01548   }
01549   return min;
01550   //cout << "complete.\n";
01551   //cout.flush();
01552 }
01553 
01554 void Container::computeGlobalEnergy (void)
01555 {
01556   // Note face intersection force is not included in energy calculation.
01557   cout << "Compute global energy..........................";
01558   cout.flush();
01559   energy=0;
01560   double dummy[3] = {0.0,0.0,0.0};
01561   // for each object in container
01562   for(o_it i=o.begin();i!=o.end();i++)
01563   {
01564     // for each vertex in object
01565     for(v_it j=i->v.begin();j!=i->v.end();j++)
01566     {
01567       // false -> do not compute force, hence dummy
01568       if(j->getFace()!=NULL)
01569       {
01570         j->setEnergy(j->getSeparationForceEnergy(dummy,false));
01571         energy+=j->getEnergy();
01572       }
01573     }
01574     // for each edge in object
01575     for(e_it j=i->e.begin();j!=i->e.end();j++)
01576     {
01577       // choice of v1 is arbitrary, v2 could have been used
01578       // false -> do not compute force, hence dummy
01579 
01580       Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
01581       j->getVertices(v1,v2,o1,o2);
01582       j->setEnergy(j->getStretchForceEnergy(v1,dummy,false)
01583             + j->getAngleForceEnergy(0,dummy,false));
01584       energy+=j->getEnergy();
01585     }
01586   }
01587   cout << "complete.\n";
01588   cout.flush();
01589 }
01590 
01591 void Container::sortAdjacentFaces (void)
01592 {
01593   // for each object
01594   for(o_it i=o.begin();i!=o.end();i++)
01595   {
01596     // for each vertex in object
01597     for(v_it j=i->v.begin();j!=i->v.end();j++)
01598     {
01599       // sort adjacent faces
01600       j->sortAdjacentFaces();
01601       //sort(j->begin(),j->end());
01602     }
01603   } 
01604 }
01605 
01606 void Container::printFrozenCP (void) const
01607 {
01608   cout << "x_coordinate y_coordinate z_coordinate state_value x_normal y_normal z_normal" << endl;
01609   for(vp_cit i=frozen.begin();i!=frozen.end();i++)
01610   {
01611     (*i)->printCP(std::cout);
01612   }
01613 }
01614 
01615 double Container::getEnergy (void) const
01616 {
01617   return energy;
01618 }
01619 
01620 void Container::addForce (double ff)
01621 {
01622   force += ff;
01623 }
01624 
01625 int Container::getFileCount (void) const
01626 {
01627   return files.size();
01628 }
01629 
01630 bool Container::vertexIsFrozen (Vertex * const v) const
01631 {
01632   return binary_search(frozen.begin(),frozen.end(),v);
01633 }
01634 
01635 void Container::updateEnergy (double d)
01636 {
01637   energy += d;
01638 }
01639 

Generated on Fri Jul 18 19:43:39 2008 for meshmorph by  doxygen 1.5.1