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
00033
00034 if (only_one == NULL)
00035 only_one = new Container();
00036
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
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
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
00103 for (o_cit i=o.begin();i!=o.end();i++)
00104 {
00105 j += i->e.size();
00106 }
00107 return j;
00108 }
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
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
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
00182 for (o_it i=o.begin();i!=o.end();i++)
00183 {
00184
00185
00186 for (v_it j=i->v.begin();j!=i->v.end();j++)
00187 {
00188
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
00195 findClosestPtToVertex(&(*j),p,sqd,ncl);
00196
00197 if(j->getFace()==ncl) continue;
00198
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
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
00282 int num_digits = i->setNumDigits();
00283
00284 map_s_ep hm;
00285
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
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
00332 str log_file = OUTPUT_DATA_DIR + NEIGHBORHOOD_FILE;
00333 std::ofstream this_file (log_file.c_str());
00334 this_file.close();
00335 }
00336
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;
00368 struct dirent *pent;
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
00375 str line = pent->d_name;
00376
00377 found = line.find(".mesh",0);
00378
00379 if (found != str::npos)
00380 {
00381
00382 files.push_back(line);
00383
00384
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
00416 for (uint count=0;count<files.size();count++)
00417 {
00418
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
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
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
00463 FILE *F = fopen(filename,"r");
00464 if (!F) { printf("Couldn't open input file %s\n",filename);return;}
00465
00466 for (char *str=fgets(line,2048,F) ; str!=NULL ; str=fgets(line,2048,F))
00467 {
00468
00469 while (strchr(" \t,",*str)!=NULL) { str++;}
00470
00471 if (strchr("V",*str)!=NULL){ v_count++; }
00472
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
00491 assessFile(obj,filename);
00492 vec_vp vp;
00493
00494 F = fopen(filename,"r");
00495 if (!F) { printf("Couldn't open input file %s\n",filename);return;}
00496
00497 for (str=fgets(line,2048,F) ; str!=NULL ; str=fgets(line,2048,F))
00498 {
00499
00500 while (strchr(" \t,",*str)!=NULL) { str++;}
00501
00502
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
00520
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
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
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
00578
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
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
00592
00593
00594 if(o2==NULL){o1 = getObjectPointer(val);}
00595 else
00596 {
00597
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
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
00624 frozen_map.insert(std::make_pair(o1,vi));
00625 }
00626 }
00627
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
00641 oi_it front;
00642
00643 s_set not_found;
00644
00645 mmap_oi frozen_map = loadFrozenMap(filename,not_found);
00646
00647
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
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
00686 for(o_cit i=o.begin();i!=o.end();i++)
00687 {
00688
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
00703 std::ofstream newfile;
00704 newfile.open (file.c_str(),std::ios::trunc);
00705 if(newfile.is_open())
00706 {
00707 newfile.precision(12);
00708
00709 for(v_cit j=i->v.begin();j!=i->v.end();j++)
00710 {
00711
00712 newfile << "Vertex "
00713 << j->getIndex() << " "
00714 << j->getCoord(0) << " "
00715 << j->getCoord(1) << " "
00716 << j->getCoord(2) << "\n";
00717 }
00718
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
00751 for (o_it i=o.begin();i!=o.end();i++)
00752 {
00753
00754 for (f_it j=i->f.begin();j!=i->f.end();j++)
00755 {
00756
00757
00758
00759
00760
00761 bp.clear();
00762 Space::instance().getBoxesToCheck(&(*j),bp);
00763
00764 if (bp.empty())
00765 {
00766 cout << "Container::checkFacesInBoxes: ERROR: NO BOXES FOR\n";
00767 j->print(std::cout);
00768 exit(1);
00769 }
00770
00771 for(bp_it k=bp.begin();k!=bp.end();k++)
00772 {
00773 bool found = false;
00774
00775 for(fp_cit m=(*k)->begin();m!=(*k)->end();m++)
00776 {
00777 if((*m)==&(*j)){found=true;break;}
00778 }
00779
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
00789
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
00806 Space::instance().clearBoxes();
00808
00809 for (o_it i=o.begin();i!=o.end();i++)
00810 {
00811
00812 for (f_it j=i->f.begin();j!=i->f.end();j++)
00813 {
00814
00815
00816
00817
00818
00819 bp.clear();
00820 Space::instance().getBoxesToCheck(&(*j),bp);
00821
00822 if (bp.empty())
00823 {
00824 cout << "ERROR: NO BOXES FOR\n";
00825 j->print(std::cout);
00826 exit(1);
00827 }
00828
00829 Space::instance().recordFace(bp,&(*j));
00830
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
00850 double a = 1.0/o.size();
00851 for (o_it i=o.begin();i!=o.end();i++)
00852 {
00853
00854 for (v_it j=i->v.begin();j!=i->v.end();j++)
00855 {
00856 if (vertexIsFrozen(&(*j))==false)
00857 {
00859
00860
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
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
00900 for (bp_it j=bp.begin();j!=bp.end();j++)
00901 {
00902
00903 for (fp_cit k=(*j)->begin();k!=(*j)->end();k++)
00904 {
00905
00906
00907
00908
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
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
00947 vec_d n = v->getNormal();
00948
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
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
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
00972
00973
00974
00975
00976
00977
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
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
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
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
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
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
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
01097
01098 vec_d sr = getSearchRegion(v);
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
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
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125 vec_bp nb;
01126 nb.reserve(VECTOR_RESERVE);
01127 Space::instance().getBoxesFromIndices(NULL,indices,nb);
01128
01129
01130
01131
01132 vec_fp cf;
01133 cf.reserve(VECTOR_RESERVE);
01134 getFacesToCheckForClosestPt(nb,v,cf);
01135
01136 sort(cf.begin(),cf.end());
01137 cf.assign(cf.begin(),unique(cf.begin(),cf.end()));
01138
01139
01140
01141
01142 Log::instance().updateClosestPtStats(cf.size(),nb.size());
01143
01144 sqd=0.0;
01145 return findClosestPtToVertexAmongFaces(v,cf,p,sqd,ncl);
01146 }
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
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
01168 if (faces.empty()==true) return gate;
01169
01170 vec_d vn = v->getNormal();
01171
01172 vec_vp verts;
01173 vec_ep edges;
01174
01175 for (fp_it j=faces.begin();j!=faces.end();j++)
01176 {
01177 assert(*j!=NULL);
01178
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
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
01203
01204 mdt_it i = tmap.begin();
01205 mdf_it j = fmap.begin();
01206 mdf_it k = j;
01207 while(i!=tmap.end())
01208 {
01209
01210
01211
01212
01213 p = (*i).second;
01214
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
01220
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
01226
01227
01228 if ( (cos_angle2 > CLOSEST_POINT_COSINE*CLOSEST_POINT_COSINE)
01229 || ((*i).first<(1.0*SCALE)))
01230 {
01231
01232 if ((*i).first<(SEARCH_RADIUS_SQ*SCALE*SCALE))
01233 {
01234
01235
01236 if ((*i).first<squareD||!squareD)
01237 {
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
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
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
01286
01287
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
01316 sort(verts.begin(),verts.end());
01317 verts.assign(verts.begin(),unique(verts.begin(),verts.end()));
01318
01319 for(vp_it i=verts.begin();i!=verts.end();i++)
01320 {
01321
01322
01323
01324
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
01330 fmap[sqd]=(*i)->getFaceNotAdjToVertex(v);
01331
01332
01333
01334
01335
01336
01337 tmap[sqd]=*((*i)->getPos());
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351 }
01352 }
01353
01354 void Container::getDistFromVertexToEdges (Vertex const * const v,vec_ep & edges,map_df & fmap,map_dt & tmap) const
01355 {
01356
01357 sort(edges.begin(),edges.end());
01358 edges.assign(edges.begin(),unique(edges.begin(),edges.end()));
01359
01360 for(ep_it i=edges.begin();i!=edges.end();i++)
01361 {
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379 vec_d p;
01380 if((*i)->getNearPtToVertex(v,p)==true)
01381 {
01382
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
01388 fmap[sqd]=(*i)->getFaceNotAdjToVertex(v);
01389
01390
01391
01392
01393
01394
01395 tmap[sqd]=p;
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420 }
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 }
01442 }
01443
01444 void Container::getDistFromVertexToFaces (Vertex const * const v,vec_fp & faces,map_df & fmap,map_dt & tmap) const
01445 {
01446
01447 sort(faces.begin(),faces.end());
01448 faces.assign(faces.begin(),unique(faces.begin(),faces.end()));
01449
01450 Point p(v->getCoord(0),v->getCoord(1),v->getCoord(2));
01451
01452 for(fp_it i=faces.begin();i!=faces.end();i++)
01453 {
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467 vec_d fn = (*i)->getNormal();
01468
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
01474 double num=dot(fn,diff);
01475
01476 bool inside = false;
01477 p.clear();
01478 if (fabs(num)>DOUBLE_EPSILON)
01479 {
01480
01481
01482
01483
01484 inside = checkNearPtOnFaceInteriorToVertex(*i,v,fn);
01485 }
01486
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
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512 fmap[p.getSqD()]=*i;
01513
01514 if(!tmap.count(p.getSqD())){tmap[p.getSqD()]=t;}
01515 }
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530 }
01531 }
01532
01533 double Container::getMinEdgeAngle (void) const
01534 {
01535
01536
01537 double min = 1E30;
01538
01539 for (o_cit i=o.begin();i!=o.end();i++)
01540 {
01541
01542 for (e_cit j=i->e.begin();j!=i->e.end();j++)
01543 {
01544
01545 double a = j->getAngle();
01546 if (a<min) min=a;
01547 }
01548 }
01549 return min;
01550
01551
01552 }
01553
01554 void Container::computeGlobalEnergy (void)
01555 {
01556
01557 cout << "Compute global energy..........................";
01558 cout.flush();
01559 energy=0;
01560 double dummy[3] = {0.0,0.0,0.0};
01561
01562 for(o_it i=o.begin();i!=o.end();i++)
01563 {
01564
01565 for(v_it j=i->v.begin();j!=i->v.end();j++)
01566 {
01567
01568 if(j->getFace()!=NULL)
01569 {
01570 j->setEnergy(j->getSeparationForceEnergy(dummy,false));
01571 energy+=j->getEnergy();
01572 }
01573 }
01574
01575 for(e_it j=i->e.begin();j!=i->e.end();j++)
01576 {
01577
01578
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
01594 for(o_it i=o.begin();i!=o.end();i++)
01595 {
01596
01597 for(v_it j=i->v.begin();j!=i->v.end();j++)
01598 {
01599
01600 j->sortAdjacentFaces();
01601
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