00001 #include "vertex.h"
00002
00003 #include <cassert>
00004 #include <cmath>
00005 #include <iostream>
00006
00007 #include "container.h"
00008 #include "edge.h"
00009 #include "face.h"
00010 #include "intersecting_faces.h"
00011 #include "misc.h"
00012 #include "nice.h"
00013 #include "object.h"
00014
00015 using std::cout;
00016 using std::endl;
00017
00018 Vertex::Vertex (const Vertex& rhs)
00019 :index(rhs.index),p(rhs.p),cl(rhs.cl),o(rhs.o),
00020 f(rhs.f),energy(rhs.energy)
00021 {
00022 }
00023
00024 Vertex& Vertex::operator = (const Vertex& rhs)
00025 {
00026 cout << "Assignment operator prohibited on instances of vertex class.\n";
00027 cout << "Vertex " << rhs.index << endl;
00028 exit(0);
00029 }
00030
00037 Vertex::Vertex (char const * triplet,Object * const q)
00038 :index(0),p(),cl(NULL),o(q),f(),energy(0)
00039 {
00040 int i;
00041 char val[80];
00042 char * eptr;
00043 char const * cp = triplet;
00044
00045
00046 while (strchr("Vertx",*triplet)!=NULL) {triplet++;}
00047
00048
00049 while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00050 i=0;
00051 while (strchr("0123456789+-eE.",*triplet)!=NULL)
00052 {
00053 val[i++] = *triplet++;
00054 }
00055 val[i]=0;
00056 index = (int) strtod(val,&eptr);
00057 if (val==eptr)
00058 {
00059 index=0;
00060 printf("Error in reading vertex index\n");
00061 return;
00062 }
00063
00064
00065 while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00066 i=0;
00067 while (strchr("0123456789+-eE.",*triplet)!=NULL)
00068 {
00069 val[i++] = *triplet++;
00070 }
00071 val[i]=0;
00072 p.push_back(strtod(val,&eptr));
00073 if (val==eptr)
00074 {
00075 printf("Error in reading vertex\n");
00076 printf("Error in reading vertex: string %s\n",cp);
00077 return;
00078 }
00079
00080
00081 while (strchr(" \t,",*triplet)!=NULL) triplet++;
00082 i=0;
00083 while (strchr("0123456789+-eE.",*triplet))
00084 {
00085 val[i++] = *triplet++;
00086 }
00087 val[i]=0;
00088 p.push_back(strtod(val,&eptr));
00089 if (val==eptr)
00090 {
00091 printf("Error in reading vertex\n");
00092 printf("Error in reading vertex: string %s\n",cp);
00093 return;
00094 }
00095
00096
00097 while (strchr(" \t,",*triplet)!=NULL) triplet++;
00098 i=0;
00099 while (strchr("0123456789+-eE.",*triplet))
00100 {
00101 val[i++] = *triplet++;
00102 }
00103 val[i]=0;
00104 p.push_back(strtod(val,&eptr));
00105 if (val==eptr)
00106 {
00107 printf("Error in reading vertex\n");
00108 printf("Error in reading vertex: string %s\n",cp);
00109 return;
00110 }
00111 }
00112
00118 int Vertex::getIndex (void) const
00119 {
00120 return index;
00121 }
00122
00128 void Vertex::print (std::ostream & target) const
00129 {
00130 char str[1024];
00131 sprintf(str,"Vertex <obj> %s <ind> %d",o->getName().c_str(),index);
00132 sprintf(str,"%s [%.15g %.15g %.15g]\n",str,p[0],p[1],p[2]);
00133 target << str;
00134 }
00135
00142 void Vertex::printCP (std::ostream &target) const
00143 {
00144 char str[1024];
00145 sprintf(str,"%.15g %.15g %.15g 1 0 0 1\n",p[0],p[1],p[2]);
00146 target << str;
00147 }
00148
00157 bool Vertex::isMatch (int const i, str const name) const
00158 {
00159 return i==index && name==o->getName();
00160 }
00168 vec_d Vertex::getNormal (void) const
00169 {
00170 double theta,thetaT=0.0,L;
00171
00172 vec_d n(3,.00);
00173
00174 for (fp_cit i=f.begin();i!=f.end();i++)
00175 {
00176
00177 vec_d t = (*i)->getNormal();
00178 L=sqrt(t[0]*t[0]+t[1]*t[1]+t[2]*t[2]);
00179
00180 theta=(*i)->getAngleProxy(this);
00181 thetaT+=theta;
00182
00183 double a = theta/L;
00184 n[0] += t[0]*a;
00185 n[1] += t[1]*a;
00186 n[2] += t[2]*a;
00187 }
00188 double a = 1.0/f.size()/thetaT;
00189 n[0] = n[0]*a;
00190 n[1] = n[1]*a;
00191 n[2] = n[2]*a;
00192 return n;
00193 }
00194
00204 double Vertex::getSeparationForceEnergy (double * const &force,bool compute_force) const
00205 {
00206 if(cl!=NULL)
00207 {
00208
00209
00210
00211 vec_fp fp(1,cl);
00212 vec_d dummy_p;
00213 double sqd=0.0;
00214 Face *dummy_ncl = NULL;
00215 Container::instance().findClosestPtToVertexAmongFaces(this,fp,dummy_p,sqd,dummy_ncl);
00216
00217
00218
00219
00220
00221
00222 vec_d s(3,0);
00223
00224
00225 double sd = sqrt(sqd);
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 double TS = TARGET_SEPARATION*SCALE;
00236
00237
00238
00239
00240 double se;
00241 if(Nice::instance().vertexIsNice(this)==false)
00242 {
00243 se=sd+TS;
00244 }
00245 else
00246 {
00247 se=sd-TS;
00248 }
00249 if(TARGET_SEPARATION==0)
00250 {
00251 cout << "Vertex::getSeparationForceEnergy: we got problems."
00252 << "TARGET_SEPARATION==0.\n";
00253 exit(0);
00254 }
00255 se = se/TARGET_SEPARATION;
00256
00257
00258
00259 if (!sd)
00260 {
00261 s = getNormal();
00262 sd = sqrt(dot(s,s));
00263 }
00264 if(compute_force==true)
00265 {
00266
00267
00268
00269 double scaled_spring_force = SEPARATION_GAIN*SEPARATION_WEIGHT*se/sd;
00270
00271
00272
00273 for (int i=0;i<3;i++) force[i]+=scaled_spring_force*s[i];
00274 }
00275
00276 return SEPARATION_WEIGHT/2.0*se*se;
00277 }
00278 else
00279 {
00280
00281
00282
00283
00284
00285
00286 double max_se = sqrt(SEARCH_RADIUS_SQ*SCALE*SCALE) + LOOP_TARGET_SEPARATION*SCALE;
00287 return SEPARATION_WEIGHT/2.0*max_se*max_se;
00288 }
00289 }
00290
00291 double Vertex::getEdgeStretchForceEnergy (double * const & force,bool flag) const
00292 {
00293 double en_ergy=0;
00294 vec_ep e;
00295 e.reserve(VECTOR_RESERVE);
00296 getAdjacentEdges(e);
00297
00298 for (ep_it j=e.begin();j!=e.end();j++)
00299 {
00300
00301 en_ergy+=(*j)->getStretchForceEnergy(this,force,flag);
00302 }
00303 return en_ergy;
00304 }
00305
00306 double Vertex::getEdgeAngleForceEnergy (double * const &force,bool flag) const
00307 {
00308
00309
00310 Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00311 vec_ep ae,oe;
00312 ae.reserve(VECTOR_RESERVE);
00313 oe.reserve(VECTOR_RESERVE);
00314
00315
00316 for(fp_cit i=f.begin();i!=f.end();i++)
00317 {
00318
00319 for(int j=0;j<3;j++)
00320 {
00321 (*i)->getEdge(j)->getVertices(v1,v2,o1,o2);
00322
00323 if(v1==this || v2==this)
00324 {
00325
00326 ae.push_back((*i)->getEdge(j));
00327 }
00328 else
00329 {
00330
00331 oe.push_back((*i)->getEdge(j));
00332 }
00333 }
00334 }
00335
00336 sort(ae.begin(),ae.end());
00337 ae.assign(ae.begin(),unique(ae.begin(),ae.end()));
00338 sort(oe.begin(),oe.end());
00339 oe.assign(oe.begin(),unique(oe.begin(),oe.end()));
00340
00341 double en_ergy=0;
00342 for (ep_it i=ae.begin();i!=ae.end();i++)
00343 {
00344
00345 en_ergy+=(*i)->getAngleReForceEnergy(force,flag);
00346 }
00347 return en_ergy;
00348 }
00349
00350 void Vertex::getAdjFaceIntForce (double * const & force) const
00351 {
00352
00353 for(fp_cit i=f.begin();i!=f.end();i++)
00354 {
00355
00356 if(Intersecting_Faces::instance().faceIsIntersectedLHS(*i)==true)
00357 {
00358
00359 Intersecting_Faces::instance().getFaceIntersectionForce(*i,force);
00360 }
00361 }
00362 }
00363
00364 void Vertex::getTotalForceEnergy (double * const &force) const
00365 {
00367
00368 getSeparationForceEnergy(force,true);
00369
00370 if (false) {cout << Container::instance().o.size() << endl;}
00371 bool print_switch = false;
00372 double fmag = 0.0;
00373 double f_old[3];
00374 f_old[0] = force[0];
00375 f_old[1] = force[1];
00376 f_old[2] = force[2];
00377 if(print_switch==true)
00378 {
00379 cout << "Vertex::getTotalForceEnergy: force ["
00380 << force[0] << " "
00381 << force[1] << " "
00382 << force[2] << "\n";
00383 }
00384
00385 if(print_switch==true)
00386 {
00387 double diff = sqrt(dot(force,force))-fmag;
00388 fmag = fmag+diff;
00389 f_old[0] = force[0];
00390 f_old[1] = force[1];
00391 f_old[2] = force[2];
00392 cout << "Vertex::getTotalForceEnergy: sep force ["
00393 << force[0] << " "
00394 << force[1] << " "
00395 << force[2] << "]"
00396 << ", diff = " << diff << endl;
00397 }
00398
00399
00401 getEdgeStretchForceEnergy(force,true);
00402
00403 if(print_switch==true)
00404 {
00405 double diff = sqrt(dot(force,force))-fmag;
00406 fmag = fmag+diff;
00407 cout << "Vertex::getTotalForceEnergy: stretch force ["
00408 << force[0]-f_old[0] << " "
00409 << force[1]-f_old[1] << " "
00410 << force[2]-f_old[2] << "]"
00411 << ", diff = " << diff << endl;
00412 f_old[0] = force[0];
00413 f_old[1] = force[1];
00414 f_old[2] = force[2];
00415 }
00416
00418 getEdgeAngleForceEnergy(force,true);
00419
00420 if(print_switch==true)
00421 {
00422 double diff = sqrt(dot(force,force))-fmag;
00423 fmag = fmag+diff;
00424 cout << "Vertex::getTotalForceEnergy: angle force ["
00425 << force[0]-f_old[0] << " "
00426 << force[1]-f_old[1] << " "
00427 << force[2]-f_old[2] << "]"
00428 << ", diff = " << diff << endl;
00429 cout << "Vertex::getTotalForceEnergy: resultant force ["
00430 << force[0] << " "
00431 << force[1] << " "
00432 << force[2] << "]\n\n";
00433 }
00434
00436 getAdjFaceIntForce(force);
00437
00438 if(print_switch==true)
00439 {
00440 double diff = sqrt(dot(force,force))-fmag;
00441 fmag = fmag+diff;
00442 cout << "Vertex::getTotalForceEnergy: face intersection force ["
00443 << force[0]-f_old[0] << " "
00444 << force[1]-f_old[1] << " "
00445 << force[2]-f_old[2] << "]"
00446 << ", diff = " << diff << endl;
00447 cout << "Vertex::getTotalForceEnergy: resultant force ["
00448 << force[0] << " "
00449 << force[1] << " "
00450 << force[2] << "]\n\n";
00451 }
00452
00453 }
00454
00455 vec_d Vertex::getNewPos (double gain)
00456 {
00457
00458 double force[3]={0,0,0};
00459 getTotalForceEnergy(force);
00460
00461 vec_d pt(3,0);
00462 for (int k=0;k<3;k++) { pt[k] = gain*force[k]+p[k]; }
00463
00464 Container::instance().addForce(sqrt(force[0]*force[0]+force[1]*force[1]+force[2]*force[2]));
00465 return pt;
00466 }
00467
00468 void Vertex::getAdjVertices (vec_vp &a) const
00469 {
00470 a.clear();
00471 vec_ep e;
00472 getAdjacentEdges(e);
00473
00474 for (ep_it i=e.begin();i!=e.end();i++)
00475 {
00476 Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00477 (*i)->getVertices(v1,v2,o1,o2);
00478
00479 if (v1!=this){a.push_back(v1);}
00480 else if (v2!=this) {a.push_back(v2);}
00481 else { printf("Error. both vertices of edge are equal to current vertex.\n"); exit(1); }
00482 }
00483 }
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 double Vertex::getSqSepDist (void) const
00499 {
00500 if (cl==NULL) return SEARCH_RADIUS_SQ;
00501
00502
00503
00504 vec_fp fp(1,cl);
00505 vec_d pp;
00506 double sqd=0.0;
00507 Face *ncl = NULL;
00508 Container::instance().findClosestPtToVertexAmongFaces(this,fp,pp,sqd,ncl);
00509
00510
00511
00512 return sqd;
00513 }
00514
00515 void Vertex::setNewPos (vec_d const & pt)
00516 {
00517 double d=0;
00519 vec_vp v;
00520 getAdjVertices(v);
00521
00522 for(vp_it i=v.begin();i!=v.end();i++)
00523 {
00524
00525 if ( !distinguishable((*i)->p[0],pt[0])&&
00526 !distinguishable((*i)->p[1],pt[1])&&
00527 !distinguishable((*i)->p[2],pt[2])){ d=DOUBLE_EPSILON*SCALE;}
00528 }
00529 p[0]=pt[0]+d;
00530 p[1]=pt[1]+d;
00531 p[2]=pt[2]+d;
00532 }
00533
00534 void Vertex::getBoundingBox (double * const & bb) const
00535 {
00536
00537 bb[0]=bb[2]=bb[4]=1E30;
00538 bb[1]=bb[3]=bb[5]=-1E30;
00539
00540 for(fp_cit i=f.begin();i!=f.end();i++)
00541 {
00542
00543 for(int j=0;j<3;j++)
00544 {
00545 Vertex *vv=(*i)->getVertex(j);
00546 if (vv->p[0]<bb[0]){bb[0]=vv->p[0];}
00547 else if (vv->p[0]>bb[1]){bb[1]=vv->p[0];}
00548 if (vv->p[1]<bb[2]){bb[2]=vv->p[1];}
00549 else if (vv->p[1]>bb[3]){bb[3]=vv->p[1];}
00550 if (vv->p[2]<bb[4]){bb[4]=vv->p[2];}
00551 else if (vv->p[2]>bb[5]){bb[5]=vv->p[2];}
00552 }
00553 }
00554 }
00555
00556 void Vertex::getAdjacentEdges (vec_ep &vec) const
00557 {
00558 vec.clear();
00559
00560 for(fp_cit i=f.begin();i!=f.end();i++)
00561 {
00562
00563 for(int j=0;j<3;j++)
00564 {
00565 Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00566 (*i)->getEdge(j)->getVertices(v1,v2,o1,o2);
00567
00568 if(v1==this || v2==this)
00569 {
00570
00571 vec.push_back((*i)->getEdge(j));
00572 }
00573 }
00574 }
00575 }
00576
00577 vec_d const * Vertex::getPos (void) const
00578 {
00579 return &p;
00580 }
00581
00582 void Vertex::setPos (double a, double b, double c)
00583 {
00584 p[0] = a;
00585 p[1] = b;
00586 p[2] = c;
00587 }
00588
00589 double Vertex::getCoord (int i) const
00590 {
00591 return p[i];
00592 }
00593
00594 Face const * Vertex::getFace (void) const
00595 {
00596 return cl;
00597 }
00598
00599 void Vertex::setFace (Face * const ff)
00600 {
00601 cl=ff;
00602 }
00603
00604 Object const * Vertex::getObject (void) const
00605 {
00606 return o;
00607 }
00608
00609 void Vertex::setObject (Object * const oo)
00610 {
00611 o=oo;
00612 }
00613
00614 fp_cit Vertex::begin (void) const
00615 {
00616 return f.begin();
00617 }
00618
00619 fp_cit Vertex::end (void) const
00620 {
00621 return f.end();
00622 }
00623
00624 void Vertex::addAdjacentFace (Face *ff)
00625 {
00626 f.push_back(ff);
00627 }
00628
00629 Face* Vertex::getAdjacentFace (int i) const
00630 {
00631 return f[i];
00632 }
00633
00634 int Vertex::getNumAdjFaces (void) const
00635 {
00636 return f.size();
00637 }
00638
00639 void Vertex::sortAdjacentFaces (void)
00640 {
00641 sort(f.begin(),f.end());
00642 }
00643
00644 double Vertex::getEnergy (void) const
00645 {
00646 return energy;
00647 }
00648
00649 void Vertex::setEnergy (double ee)
00650 {
00651 energy = ee;
00652 }
00653
00654 bool Vertex::faceIsAdjacent (Face const * const ff) const
00655 {
00656 return binary_search(f.begin(),f.end(),ff);
00657 }
00658
00659 double Vertex::getSqVirtualDisp (double gain)
00660 {
00661
00662
00663 vec_d q = getNewPos(gain);
00664 return (q[0]-getCoord(0))*(q[0]-getCoord(0))+
00665 (q[1]-getCoord(1))*(q[1]-getCoord(1))+
00666 (q[2]-getCoord(2))*(q[2]-getCoord(2));
00667 }
00668
00669 Face* Vertex::getFaceNotAdjToVertex (Vertex const * const cv) const
00670 {
00671
00672 for(fp_cit i=begin();i!=end();i++)
00673 {
00674
00675 if(binary_search(cv->begin(),cv->end(),*i)==false)
00676 {
00677 return *i;
00678 }
00679 }
00680 cout << "Vertex::getFaceNotAdjToVertex: Error:"
00681 << "no acceptable face was found!\n";
00682 exit(0);
00683 }
00684