Changeset 3d919e
- Timestamp:
- Jul 7, 2009, 6:49:32 AM (15 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- 260b2f
- Parents:
- ca2587
- Location:
- src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundary.cpp
rca2587 r3d919e 1 1 #include "boundary.hpp" 2 #include "linkedcell.hpp" 3 #include "molecules.hpp" 4 #include <gsl/gsl_matrix.h> 5 #include <gsl/gsl_linalg.h> 6 #include <gsl/gsl_multimin.h> 7 #include <gsl/gsl_permutation.h> 2 8 3 9 #define DEBUG 1 4 #define DoSingleStepOutput 010 #define DoSingleStepOutput 1 5 11 #define DoTecplotOutput 1 6 12 #define DoRaster3DOutput 1 … … 9 15 #define Raster3DSuffix ".r3d" 10 16 #define VRMLSUffix ".wrl" 11 #define HULLEPSILON MYEPSILON17 #define HULLEPSILON 1e-7 12 18 13 19 // ======================================== Points on Boundary ================================= … … 15 21 BoundaryPointSet::BoundaryPointSet() 16 22 { 17 18 23 LinesCount = 0; 24 Nr = -1; 19 25 } 20 26 ; … … 22 28 BoundaryPointSet::BoundaryPointSet(atom *Walker) 23 29 { 24 25 26 30 node = Walker; 31 LinesCount = 0; 32 Nr = Walker->nr; 27 33 } 28 34 ; … … 30 36 BoundaryPointSet::~BoundaryPointSet() 31 37 { 32 33 34 35 38 cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl; 39 if (!lines.empty()) 40 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl; 41 node = NULL; 36 42 } 37 43 ; … … 39 45 void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 40 46 { 41 42 43 44 45 46 47 48 49 50 51 47 cout << Verbose(6) << "Adding " << *this << " to line " << *line << "." 48 << endl; 49 if (line->endpoints[0] == this) 50 { 51 lines.insert(LinePair(line->endpoints[1]->Nr, line)); 52 } 53 else 54 { 55 lines.insert(LinePair(line->endpoints[0]->Nr, line)); 56 } 57 LinesCount++; 52 58 } 53 59 ; … … 56 62 operator <<(ostream &ost, BoundaryPointSet &a) 57 63 { 58 59 64 ost << "[" << a.Nr << "|" << a.node->Name << "]"; 65 return ost; 60 66 } 61 67 ; … … 65 71 BoundaryLineSet::BoundaryLineSet() 66 72 { 67 68 69 70 73 for (int i = 0; i < 2; i++) 74 endpoints[i] = NULL; 75 TrianglesCount = 0; 76 Nr = -1; 71 77 } 72 78 ; … … 74 80 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number) 75 81 { 76 77 78 79 80 81 82 83 84 85 82 // set number 83 Nr = number; 84 // set endpoints in ascending order 85 SetEndpointsOrdered(endpoints, Point[0], Point[1]); 86 // add this line to the hash maps of both endpoints 87 Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding. 88 Point[1]->AddLine(this); // 89 // clear triangles list 90 TrianglesCount = 0; 91 cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl; 86 92 } 87 93 ; … … 89 95 BoundaryLineSet::~BoundaryLineSet() 90 96 { 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 97 int Numbers[2]; 98 Numbers[0] = endpoints[1]->Nr; 99 Numbers[1] = endpoints[0]->Nr; 100 for (int i = 0; i < 2; i++) { 101 cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 102 endpoints[i]->lines.erase(Numbers[i]); 103 if (endpoints[i]->lines.empty()) { 104 cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl; 105 if (endpoints[i] != NULL) { 106 delete(endpoints[i]); 107 endpoints[i] = NULL; 108 } else 109 cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl; 110 } else 111 cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl; 112 } 113 if (!triangles.empty()) 114 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl; 109 115 } 110 116 ; … … 113 119 BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 114 120 { 115 116 117 118 121 cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." 122 << endl; 123 triangles.insert(TrianglePair(triangle->Nr, triangle)); 124 TrianglesCount++; 119 125 } 120 126 ; … … 123 129 operator <<(ostream &ost, BoundaryLineSet &a) 124 130 { 125 126 127 131 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," 132 << a.endpoints[1]->node->Name << "]"; 133 return ost; 128 134 } 129 135 ; … … 134 140 BoundaryTriangleSet::BoundaryTriangleSet() 135 141 { 136 137 138 139 140 141 142 for (int i = 0; i < 3; i++) 143 { 144 endpoints[i] = NULL; 145 lines[i] = NULL; 146 } 147 Nr = -1; 142 148 } 143 149 ; 144 150 145 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], 146 int number) 147 { 148 // set number 149 Nr = number; 150 // set lines 151 cout << Verbose(5) << "New triangle " << Nr << ":" << endl; 152 for (int i = 0; i < 3; i++) 153 { 154 lines[i] = line[i]; 155 lines[i]->AddTriangle(this); 156 } 157 // get ascending order of endpoints 158 map<int, class BoundaryPointSet *> OrderMap; 159 for (int i = 0; i < 3; i++) 160 // for all three lines 161 for (int j = 0; j < 2; j++) 162 { // for both endpoints 163 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 164 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 165 // and we don't care whether insertion fails 166 } 167 // set endpoints 168 int Counter = 0; 169 cout << Verbose(6) << " with end points "; 170 for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner 171 != OrderMap.end(); runner++) 172 { 173 endpoints[Counter] = runner->second; 174 cout << " " << *endpoints[Counter]; 175 Counter++; 176 } 177 if (Counter < 3) 178 { 179 cerr << "ERROR! We have a triangle with only two distinct endpoints!" 180 << endl; 181 //exit(1); 182 } 183 cout << "." << endl; 151 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) 152 { 153 // set number 154 Nr = number; 155 // set lines 156 cout << Verbose(5) << "New triangle " << Nr << ":" << endl; 157 for (int i = 0; i < 3; i++) 158 { 159 lines[i] = line[i]; 160 lines[i]->AddTriangle(this); 161 } 162 // get ascending order of endpoints 163 map<int, class BoundaryPointSet *> OrderMap; 164 for (int i = 0; i < 3; i++) 165 // for all three lines 166 for (int j = 0; j < 2; j++) 167 { // for both endpoints 168 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 169 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 170 // and we don't care whether insertion fails 171 } 172 // set endpoints 173 int Counter = 0; 174 cout << Verbose(6) << " with end points "; 175 for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner 176 != OrderMap.end(); runner++) 177 { 178 endpoints[Counter] = runner->second; 179 cout << " " << *endpoints[Counter]; 180 Counter++; 181 } 182 if (Counter < 3) 183 { 184 cerr << "ERROR! We have a triangle with only two distinct endpoints!" 185 << endl; 186 //exit(1); 187 } 188 cout << "." << endl; 184 189 } 185 190 ; … … 187 192 BoundaryTriangleSet::~BoundaryTriangleSet() 188 193 { 189 190 191 192 193 194 195 196 197 198 199 200 201 194 for (int i = 0; i < 3; i++) { 195 cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl; 196 lines[i]->triangles.erase(Nr); 197 if (lines[i]->triangles.empty()) { 198 cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl; 199 if (lines[i] != NULL) { 200 delete (lines[i]); 201 lines[i] = NULL; 202 } else 203 cerr << "ERROR: This line " << i << " has already been free'd." << endl; 204 } else 205 cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl; 206 } 202 207 } 203 208 ; … … 206 211 BoundaryTriangleSet::GetNormalVector(Vector &OtherVector) 207 212 { 208 209 210 211 212 213 214 213 // get normal vector 214 NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x, 215 &endpoints[2]->node->x); 216 217 // make it always point inward (any offset vector onto plane projected onto normal vector suffices) 218 if (NormalVector.Projection(&OtherVector) > 0) 219 NormalVector.Scale(-1.); 215 220 } 216 221 ; … … 219 224 operator <<(ostream &ost, BoundaryTriangleSet &a) 220 225 { 221 222 223 226 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," 227 << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]"; 228 return ost; 224 229 } 225 230 ; 231 232 233 // ============================ CandidateForTesselation ============================= 234 235 CandidateForTesselation::CandidateForTesselation( 236 atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter 237 ) { 238 point = candidate; 239 BaseLine = line; 240 OptCenter.CopyVector(&OptCandidateCenter); 241 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); 242 } 243 244 CandidateForTesselation::~CandidateForTesselation() { 245 point = NULL; 246 } 226 247 227 248 // ========================================== F U N C T I O N S ================================= … … 235 256 GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2) 236 257 { 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 258 class BoundaryLineSet * lines[2] = 259 { line1, line2 }; 260 class BoundaryPointSet *node = NULL; 261 map<int, class BoundaryPointSet *> OrderMap; 262 pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest; 263 for (int i = 0; i < 2; i++) 264 // for both lines 265 for (int j = 0; j < 2; j++) 266 { // for both endpoints 267 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> ( 268 lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j])); 269 if (!OrderTest.second) 270 { // if insertion fails, we have common endpoint 271 node = OrderTest.first->second; 272 cout << Verbose(5) << "Common endpoint of lines " << *line1 273 << " and " << *line2 << " is: " << *node << "." << endl; 274 j = 2; 275 i = 2; 276 break; 277 } 278 } 279 return node; 259 280 } 260 281 ; … … 270 291 GetBoundaryPoints(ofstream *out, molecule *mol) 271 292 { 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 //*out << Verbose(1) << "Axisvector is ";292 //AxisVector.Output(out);293 //*out << " and AngleReferenceVector is ";294 //AngleReferenceVector.Output(out);295 //*out << "." << endl;296 //*out << " and AngleReferenceNormalVector is ";297 //AngleReferenceNormalVector.Output(out);298 //*out << "." << endl;299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 //{369 //*out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;370 //int i=0;371 //for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {372 //if (runner != BoundaryPoints[axis].begin())373 //*out << ", " << i << ": " << *runner->second.second;374 //else375 //*out << i << ": " << *runner->second.second;376 //i++;377 //}378 //*out << endl;379 //}380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 //*out << "SideA: ";417 //SideA.Output(out);418 //*out << endl;419 420 421 422 //*out << "SideB: ";423 //SideB.Output(out);424 //*out << endl;425 426 427 428 429 //*out << "SideC: ";430 //SideC.Output(out);431 //*out << endl;432 433 434 435 //*out << "SideH: ";436 //SideH.Output(out);437 //*out << endl;438 439 440 441 442 443 444 445 446 447 448 449 450 451 //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 293 atom *Walker = NULL; 294 PointMap PointsOnBoundary; 295 LineMap LinesOnBoundary; 296 TriangleMap TrianglesOnBoundary; 297 298 *out << Verbose(1) << "Finding all boundary points." << endl; 299 Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr) 300 BoundariesTestPair BoundaryTestPair; 301 Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector; 302 double radius, angle; 303 // 3a. Go through every axis 304 for (int axis = 0; axis < NDIM; axis++) 305 { 306 AxisVector.Zero(); 307 AngleReferenceVector.Zero(); 308 AngleReferenceNormalVector.Zero(); 309 AxisVector.x[axis] = 1.; 310 AngleReferenceVector.x[(axis + 1) % NDIM] = 1.; 311 AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.; 312 // *out << Verbose(1) << "Axisvector is "; 313 // AxisVector.Output(out); 314 // *out << " and AngleReferenceVector is "; 315 // AngleReferenceVector.Output(out); 316 // *out << "." << endl; 317 // *out << " and AngleReferenceNormalVector is "; 318 // AngleReferenceNormalVector.Output(out); 319 // *out << "." << endl; 320 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours 321 Walker = mol->start; 322 while (Walker->next != mol->end) 323 { 324 Walker = Walker->next; 325 Vector ProjectedVector; 326 ProjectedVector.CopyVector(&Walker->x); 327 ProjectedVector.ProjectOntoPlane(&AxisVector); 328 // correct for negative side 329 //if (Projection(y) < 0) 330 //angle = 2.*M_PI - angle; 331 radius = ProjectedVector.Norm(); 332 if (fabs(radius) > MYEPSILON) 333 angle = ProjectedVector.Angle(&AngleReferenceVector); 334 else 335 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues 336 337 //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl; 338 if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) 339 { 340 angle = 2. * M_PI - angle; 341 } 342 //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): "; 343 //ProjectedVector.Output(out); 344 //*out << endl; 345 BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, 346 DistancePair (radius, Walker))); 347 if (BoundaryTestPair.second) 348 { // successfully inserted 349 } 350 else 351 { // same point exists, check first r, then distance of original vectors to center of gravity 352 *out << Verbose(2) 353 << "Encountered two vectors whose projection onto axis " 354 << axis << " is equal: " << endl; 355 *out << Verbose(2) << "Present vector: "; 356 BoundaryTestPair.first->second.second->x.Output(out); 357 *out << endl; 358 *out << Verbose(2) << "New vector: "; 359 Walker->x.Output(out); 360 *out << endl; 361 double tmp = ProjectedVector.Norm(); 362 if (tmp > BoundaryTestPair.first->second.first) 363 { 364 BoundaryTestPair.first->second.first = tmp; 365 BoundaryTestPair.first->second.second = Walker; 366 *out << Verbose(2) << "Keeping new vector." << endl; 367 } 368 else if (tmp == BoundaryTestPair.first->second.first) 369 { 370 if (BoundaryTestPair.first->second.second->x.ScalarProduct( 371 &BoundaryTestPair.first->second.second->x) 372 < Walker->x.ScalarProduct(&Walker->x)) 373 { // Norm() does a sqrt, which makes it a lot slower 374 BoundaryTestPair.first->second.second = Walker; 375 *out << Verbose(2) << "Keeping new vector." << endl; 376 } 377 else 378 { 379 *out << Verbose(2) << "Keeping present vector." << endl; 380 } 381 } 382 else 383 { 384 *out << Verbose(2) << "Keeping present vector." << endl; 385 } 386 } 387 } 388 // printing all inserted for debugging 389 // { 390 // *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl; 391 // int i=0; 392 // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) { 393 // if (runner != BoundaryPoints[axis].begin()) 394 // *out << ", " << i << ": " << *runner->second.second; 395 // else 396 // *out << i << ": " << *runner->second.second; 397 // i++; 398 // } 399 // *out << endl; 400 // } 401 // 3c. throw out points whose distance is less than the mean of left and right neighbours 402 bool flag = false; 403 do 404 { // do as long as we still throw one out per round 405 *out << Verbose(1) 406 << "Looking for candidates to kick out by convex condition ... " 407 << endl; 408 flag = false; 409 Boundaries::iterator left = BoundaryPoints[axis].end(); 410 Boundaries::iterator right = BoundaryPoints[axis].end(); 411 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner 412 != BoundaryPoints[axis].end(); runner++) 413 { 414 // set neighbours correctly 415 if (runner == BoundaryPoints[axis].begin()) 416 { 417 left = BoundaryPoints[axis].end(); 418 } 419 else 420 { 421 left = runner; 422 } 423 left--; 424 right = runner; 425 right++; 426 if (right == BoundaryPoints[axis].end()) 427 { 428 right = BoundaryPoints[axis].begin(); 429 } 430 // check distance 431 432 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector) 433 { 434 Vector SideA, SideB, SideC, SideH; 435 SideA.CopyVector(&left->second.second->x); 436 SideA.ProjectOntoPlane(&AxisVector); 437 // *out << "SideA: "; 438 // SideA.Output(out); 439 // *out << endl; 440 441 SideB.CopyVector(&right->second.second->x); 442 SideB.ProjectOntoPlane(&AxisVector); 443 // *out << "SideB: "; 444 // SideB.Output(out); 445 // *out << endl; 446 447 SideC.CopyVector(&left->second.second->x); 448 SideC.SubtractVector(&right->second.second->x); 449 SideC.ProjectOntoPlane(&AxisVector); 450 // *out << "SideC: "; 451 // SideC.Output(out); 452 // *out << endl; 453 454 SideH.CopyVector(&runner->second.second->x); 455 SideH.ProjectOntoPlane(&AxisVector); 456 // *out << "SideH: "; 457 // SideH.Output(out); 458 // *out << endl; 459 460 // calculate each length 461 double a = SideA.Norm(); 462 //double b = SideB.Norm(); 463 //double c = SideC.Norm(); 464 double h = SideH.Norm(); 465 // calculate the angles 466 double alpha = SideA.Angle(&SideH); 467 double beta = SideA.Angle(&SideC); 468 double gamma = SideB.Angle(&SideH); 469 double delta = SideC.Angle(&SideH); 470 double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha 471 < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.); 472 // *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl; 473 //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl; 474 if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) 475 < MYEPSILON) && (h < MinDistance)) 476 { 477 // throw out point 478 //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl; 479 BoundaryPoints[axis].erase(runner); 480 flag = true; 481 } 482 } 483 } 484 } 485 while (flag); 486 } 487 return BoundaryPoints; 467 488 } 468 489 ; … … 478 499 double * 479 500 GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, 480 481 { 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 501 bool IsAngstroem) 502 { 503 // get points on boundary of NULL was given as parameter 504 bool BoundaryFreeFlag = false; 505 Boundaries *BoundaryPoints = BoundaryPtr; 506 if (BoundaryPoints == NULL) 507 { 508 BoundaryFreeFlag = true; 509 BoundaryPoints = GetBoundaryPoints(out, mol); 510 } 511 else 512 { 513 *out << Verbose(1) << "Using given boundary points set." << endl; 514 } 515 // determine biggest "diameter" of cluster for each axis 516 Boundaries::iterator Neighbour, OtherNeighbour; 517 double *GreatestDiameter = new double[NDIM]; 518 for (int i = 0; i < NDIM; i++) 519 GreatestDiameter[i] = 0.; 520 double OldComponent, tmp, w1, w2; 521 Vector DistanceVector, OtherVector; 522 int component, Othercomponent; 523 for (int axis = 0; axis < NDIM; axis++) 524 { // regard each projected plane 525 //*out << Verbose(1) << "Current axis is " << axis << "." << endl; 526 for (int j = 0; j < 2; j++) 527 { // and for both axis on the current plane 528 component = (axis + j + 1) % NDIM; 529 Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM; 530 //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl; 531 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner 532 != BoundaryPoints[axis].end(); runner++) 533 { 534 //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl; 535 // seek for the neighbours pair where the Othercomponent sign flips 536 Neighbour = runner; 537 Neighbour++; 538 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around 539 Neighbour = BoundaryPoints[axis].begin(); 540 DistanceVector.CopyVector(&runner->second.second->x); 541 DistanceVector.SubtractVector(&Neighbour->second.second->x); 542 do 543 { // seek for neighbour pair where it flips 544 OldComponent = DistanceVector.x[Othercomponent]; 545 Neighbour++; 546 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around 547 Neighbour = BoundaryPoints[axis].begin(); 548 DistanceVector.CopyVector(&runner->second.second->x); 549 DistanceVector.SubtractVector(&Neighbour->second.second->x); 550 //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl; 551 } 552 while ((runner != Neighbour) && (fabs(OldComponent / fabs( 553 OldComponent) - DistanceVector.x[Othercomponent] / fabs( 554 DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip 555 if (runner != Neighbour) 556 { 557 OtherNeighbour = Neighbour; 558 if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around 559 OtherNeighbour = BoundaryPoints[axis].end(); 560 OtherNeighbour--; 561 //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 562 // now we have found the pair: Neighbour and OtherNeighbour 563 OtherVector.CopyVector(&runner->second.second->x); 564 OtherVector.SubtractVector(&OtherNeighbour->second.second->x); 565 //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl; 566 //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl; 567 // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour 568 w1 = fabs(OtherVector.x[Othercomponent]); 569 w2 = fabs(DistanceVector.x[Othercomponent]); 570 tmp = fabs((w1 * DistanceVector.x[component] + w2 571 * OtherVector.x[component]) / (w1 + w2)); 572 // mark if it has greater diameter 573 //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl; 574 GreatestDiameter[component] = (GreatestDiameter[component] 575 > tmp) ? GreatestDiameter[component] : tmp; 576 } //else 577 //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl; 578 } 579 } 580 } 581 *out << Verbose(0) << "RESULT: The biggest diameters are " 582 << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " 583 << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" 584 : "atomiclength") << "." << endl; 585 586 // free reference lists 587 if (BoundaryFreeFlag) 588 delete[] (BoundaryPoints); 589 590 return GreatestDiameter; 570 591 } 571 592 ; … … 579 600 void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol) 580 601 { 581 582 583 584 585 586 587 588 589 590 591 592 *vrmlfile << "Sphere {" << endl << ""; // 2 is sphere type593 594 595 596 597 598 599 600 601 *vrmlfile << "3" << endl << ""; // 2 is round-ended cylinder type602 603 604 605 606 607 608 609 610 611 612 *vrmlfile << "1" << endl << ""; // 1 is triangle type613 614 for (int j=0;j<NDIM;j++)// and for each node all NDIM coordinates615 616 617 618 *vrmlfile << "1. 0. 0." << endl;// red as colour619 *vrmlfile << "18" << endl << "0.5 0.5 0.5" << endl; // 18 is transparency type for previous object620 621 622 623 624 602 atom *Walker = mol->start; 603 bond *Binder = mol->first; 604 int i; 605 Vector *center = mol->DetermineCenterOfAll(out); 606 if (vrmlfile != NULL) { 607 //cout << Verbose(1) << "Writing Raster3D file ... "; 608 *vrmlfile << "#VRML V2.0 utf8" << endl; 609 *vrmlfile << "#Created by molecuilder" << endl; 610 *vrmlfile << "#All atoms as spheres" << endl; 611 while (Walker->next != mol->end) { 612 Walker = Walker->next; 613 *vrmlfile << "Sphere {" << endl << " "; // 2 is sphere type 614 for (i=0;i<NDIM;i++) 615 *vrmlfile << Walker->x.x[i]+center->x[i] << " "; 616 *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour 617 } 618 619 *vrmlfile << "# All bonds as vertices" << endl; 620 while (Binder->next != mol->last) { 621 Binder = Binder->next; 622 *vrmlfile << "3" << endl << " "; // 2 is round-ended cylinder type 623 for (i=0;i<NDIM;i++) 624 *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " "; 625 *vrmlfile << "\t0.03\t"; 626 for (i=0;i<NDIM;i++) 627 *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " "; 628 *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour 629 } 630 631 *vrmlfile << "# All tesselation triangles" << endl; 632 for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) { 633 *vrmlfile << "1" << endl << " "; // 1 is triangle type 634 for (i=0;i<3;i++) { // print each node 635 for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates 636 *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " "; 637 *vrmlfile << "\t"; 638 } 639 *vrmlfile << "1. 0. 0." << endl; // red as colour 640 *vrmlfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object 641 } 642 } else { 643 cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl; 644 } 645 delete(center); 625 646 }; 626 647 … … 633 654 void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol) 634 655 { 635 636 637 638 639 640 641 642 643 644 645 646 *rasterfile << "2" << endl << " ";// 2 is sphere type647 648 649 650 651 652 653 654 655 *rasterfile << "3" << endl << " ";// 2 is round-ended cylinder type656 657 658 659 660 661 662 663 664 665 *rasterfile << "8\n 25. -1. 1. 1. 1. 0.0 0 0 0 2\n SOLID 1.0 0.0 0.0\n BACKFACE 0.3 0.3 1.00 0\n";666 667 *rasterfile << "1" << endl << " ";// 1 is triangle type668 for (i=0;i<3;i++) {// print each node669 for (int j=0;j<NDIM;j++)// and for each node all NDIM coordinates670 671 672 673 *rasterfile << "1. 0. 0." << endl;// red as colour674 //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl;// 18 is transparency type for previous object675 676 *rasterfile << "9\nterminating special property\n";677 678 679 680 656 atom *Walker = mol->start; 657 bond *Binder = mol->first; 658 int i; 659 Vector *center = mol->DetermineCenterOfAll(out); 660 if (rasterfile != NULL) { 661 //cout << Verbose(1) << "Writing Raster3D file ... "; 662 *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl; 663 *rasterfile << "@header.r3d" << endl; 664 *rasterfile << "# All atoms as spheres" << endl; 665 while (Walker->next != mol->end) { 666 Walker = Walker->next; 667 *rasterfile << "2" << endl << " "; // 2 is sphere type 668 for (i=0;i<NDIM;i++) 669 *rasterfile << Walker->x.x[i]+center->x[i] << " "; 670 *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour 671 } 672 673 *rasterfile << "# All bonds as vertices" << endl; 674 while (Binder->next != mol->last) { 675 Binder = Binder->next; 676 *rasterfile << "3" << endl << " "; // 2 is round-ended cylinder type 677 for (i=0;i<NDIM;i++) 678 *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " "; 679 *rasterfile << "\t0.03\t"; 680 for (i=0;i<NDIM;i++) 681 *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " "; 682 *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour 683 } 684 685 *rasterfile << "# All tesselation triangles" << endl; 686 *rasterfile << "8\n 25. -1. 1. 1. 1. 0.0 0 0 0 2\n SOLID 1.0 0.0 0.0\n BACKFACE 0.3 0.3 1.0 0 0\n"; 687 for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) { 688 *rasterfile << "1" << endl << " "; // 1 is triangle type 689 for (i=0;i<3;i++) { // print each node 690 for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates 691 *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " "; 692 *rasterfile << "\t"; 693 } 694 *rasterfile << "1. 0. 0." << endl; // red as colour 695 //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object 696 } 697 *rasterfile << "9\n terminating special property\n"; 698 } else { 699 cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl; 700 } 701 delete(center); 681 702 }; 682 703 … … 688 709 void 689 710 write_tecplot_file(ofstream *out, ofstream *tecplot, 690 691 { 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 711 class Tesselation *TesselStruct, class molecule *mol, int N) 712 { 713 if (tecplot != NULL) 714 { 715 *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl; 716 *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl; 717 *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" 718 << TesselStruct->PointsOnBoundaryCount << ", E=" 719 << TesselStruct->TrianglesOnBoundaryCount 720 << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl; 721 int *LookupList = new int[mol->AtomCount]; 722 for (int i = 0; i < mol->AtomCount; i++) 723 LookupList[i] = -1; 724 725 // print atom coordinates 726 *out << Verbose(2) << "The following triangles were created:"; 727 int Counter = 1; 728 atom *Walker = NULL; 729 for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target 730 != TesselStruct->PointsOnBoundary.end(); target++) 731 { 732 Walker = target->second->node; 733 LookupList[Walker->nr] = Counter++; 734 *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " 735 << Walker->x.x[2] << " " << endl; 736 } 737 *tecplot << endl; 738 // print connectivity 739 for (TriangleMap::iterator runner = 740 TesselStruct->TrianglesOnBoundary.begin(); runner 741 != TesselStruct->TrianglesOnBoundary.end(); runner++) 742 { 743 *out << " " << runner->second->endpoints[0]->node->Name << "<->" 744 << runner->second->endpoints[1]->node->Name << "<->" 745 << runner->second->endpoints[2]->node->Name; 746 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " 747 << LookupList[runner->second->endpoints[1]->node->nr] << " " 748 << LookupList[runner->second->endpoints[2]->node->nr] << endl; 749 } 750 delete[] (LookupList); 751 *out << endl; 752 } 732 753 } 733 754 … … 743 764 double 744 765 VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration, 745 746 { 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 //*out << Verbose(1) << "Listing PointsOnBoundary:";798 //for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {799 //*out << " " << *runner->second;800 //}801 //*out << endl;802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 766 Boundaries *BoundaryPtr, molecule *mol) 767 { 768 bool IsAngstroem = configuration->GetIsAngstroem(); 769 atom *Walker = NULL; 770 struct Tesselation *TesselStruct = new Tesselation; 771 bool BoundaryFreeFlag = false; 772 Boundaries *BoundaryPoints = BoundaryPtr; 773 double volume = 0.; 774 double PyramidVolume = 0.; 775 double G, h; 776 Vector x, y; 777 double a, b, c; 778 779 //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line. 780 781 // 1. calculate center of gravity 782 *out << endl; 783 Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out); 784 785 // 2. translate all points into CoG 786 *out << Verbose(1) << "Translating system to Center of Gravity." << endl; 787 Walker = mol->start; 788 while (Walker->next != mol->end) 789 { 790 Walker = Walker->next; 791 Walker->x.Translate(CenterOfGravity); 792 } 793 794 // 3. Find all points on the boundary 795 if (BoundaryPoints == NULL) 796 { 797 BoundaryFreeFlag = true; 798 BoundaryPoints = GetBoundaryPoints(out, mol); 799 } 800 else 801 { 802 *out << Verbose(1) << "Using given boundary points set." << endl; 803 } 804 805 // 4. fill the boundary point list 806 for (int axis = 0; axis < NDIM; axis++) 807 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner 808 != BoundaryPoints[axis].end(); runner++) 809 { 810 TesselStruct->AddPoint(runner->second.second); 811 } 812 813 *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount 814 << " points on the convex boundary." << endl; 815 // now we have the whole set of edge points in the BoundaryList 816 817 // listing for debugging 818 // *out << Verbose(1) << "Listing PointsOnBoundary:"; 819 // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) { 820 // *out << " " << *runner->second; 821 // } 822 // *out << endl; 823 824 // 5a. guess starting triangle 825 TesselStruct->GuessStartingTriangle(out); 826 827 // 5b. go through all lines, that are not yet part of two triangles (only of one so far) 828 TesselStruct->TesselateOnBoundary(out, configuration, mol); 829 830 *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount 831 << " triangles with " << TesselStruct->LinesOnBoundaryCount 832 << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." 833 << endl; 834 835 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes 836 *out << Verbose(1) 837 << "Calculating the volume of the pyramids formed out of triangles and center of gravity." 838 << endl; 839 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner 840 != TesselStruct->TrianglesOnBoundary.end(); runner++) 841 { // go through every triangle, calculate volume of its pyramid with CoG as peak 842 x.CopyVector(&runner->second->endpoints[0]->node->x); 843 x.SubtractVector(&runner->second->endpoints[1]->node->x); 844 y.CopyVector(&runner->second->endpoints[0]->node->x); 845 y.SubtractVector(&runner->second->endpoints[2]->node->x); 846 a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared( 847 &runner->second->endpoints[1]->node->x)); 848 b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared( 849 &runner->second->endpoints[2]->node->x)); 850 c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared( 851 &runner->second->endpoints[1]->node->x)); 852 G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle 853 x.MakeNormalVector(&runner->second->endpoints[0]->node->x, 854 &runner->second->endpoints[1]->node->x, 855 &runner->second->endpoints[2]->node->x); 856 x.Scale(runner->second->endpoints[1]->node->x.Projection(&x)); 857 h = x.Norm(); // distance of CoG to triangle 858 PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak) 859 *out << Verbose(2) << "Area of triangle is " << G << " " 860 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " 861 << h << " and the volume is " << PyramidVolume << " " 862 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl; 863 volume += PyramidVolume; 864 } 865 *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) 866 << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." 867 << endl; 868 869 // 7. translate all points back from CoG 870 *out << Verbose(1) << "Translating system back from Center of Gravity." 871 << endl; 872 CenterOfGravity->Scale(-1); 873 Walker = mol->start; 874 while (Walker->next != mol->end) 875 { 876 Walker = Walker->next; 877 Walker->x.Translate(CenterOfGravity); 878 } 879 880 // 8. Store triangles in tecplot file 881 string OutputName(filename); 882 OutputName.append(TecplotSuffix); 883 ofstream *tecplot = new ofstream(OutputName.c_str()); 884 write_tecplot_file(out, tecplot, TesselStruct, mol, 0); 885 tecplot->close(); 886 delete(tecplot); 887 888 // free reference lists 889 if (BoundaryFreeFlag) 890 delete[] (BoundaryPoints); 891 892 return volume; 872 893 } 873 894 ; … … 883 904 void 884 905 PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, 885 886 { 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 906 double ClusterVolume, double celldensity) 907 { 908 // transform to PAS 909 mol->PrincipalAxisSystem(out, true); 910 911 // some preparations beforehand 912 bool IsAngstroem = configuration->GetIsAngstroem(); 913 Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol); 914 double clustervolume; 915 if (ClusterVolume == 0) 916 clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration, 917 BoundaryPoints, mol); 918 else 919 clustervolume = ClusterVolume; 920 double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, 921 IsAngstroem); 922 Vector BoxLengths; 923 int repetition[NDIM] = 924 { 1, 1, 1 }; 925 int TotalNoClusters = 1; 926 for (int i = 0; i < NDIM; i++) 927 TotalNoClusters *= repetition[i]; 928 929 // sum up the atomic masses 930 double totalmass = 0.; 931 atom *Walker = mol->start; 932 while (Walker->next != mol->end) 933 { 934 Walker = Walker->next; 935 totalmass += Walker->type->mass; 936 } 937 *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) 938 << totalmass << " atomicmassunit." << endl; 939 940 *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) 941 << totalmass / clustervolume << " atomicmassunit/" 942 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl; 943 944 // solve cubic polynomial 945 *out << Verbose(1) << "Solving equidistant suspension in water problem ..." 946 << endl; 947 double cellvolume; 948 if (IsAngstroem) 949 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass 950 / clustervolume)) / (celldensity - 1); 951 else 952 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass 953 / clustervolume)) / (celldensity - 1); 954 *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity 955 << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" 956 : "atomiclength") << "^3." << endl; 957 958 double minimumvolume = TotalNoClusters * (GreatestDiameter[0] 959 * GreatestDiameter[1] * GreatestDiameter[2]); 960 *out << Verbose(1) 961 << "Minimum volume of the convex envelope contained in a rectangular box is " 962 << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" 963 : "atomiclength") << "^3." << endl; 964 if (minimumvolume > cellvolume) 965 { 966 cerr << Verbose(0) 967 << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" 968 << endl; 969 cout << Verbose(0) 970 << "Setting Box dimensions to minimum possible, the greatest diameters." 971 << endl; 972 for (int i = 0; i < NDIM; i++) 973 BoxLengths.x[i] = GreatestDiameter[i]; 974 mol->CenterEdge(out, &BoxLengths); 975 } 976 else 977 { 978 BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] 979 * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]); 980 BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] 981 * GreatestDiameter[1] + repetition[0] * repetition[2] 982 * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] 983 * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]); 984 BoxLengths.x[2] = minimumvolume - cellvolume; 985 double x0 = 0., x1 = 0., x2 = 0.; 986 if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], 987 BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return 988 *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 989 << " ." << endl; 990 else 991 { 992 *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 993 << " and " << x1 << " and " << x2 << " ." << endl; 994 x0 = x2; // sorted in ascending order 995 } 996 997 cellvolume = 1; 998 for (int i = 0; i < NDIM; i++) 999 { 1000 BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]); 1001 cellvolume *= BoxLengths.x[i]; 1002 } 1003 1004 // set new box dimensions 1005 *out << Verbose(0) << "Translating to box with these boundaries." << endl; 1006 mol->CenterInBox((ofstream *) &cout, &BoxLengths); 1007 } 1008 // update Box of atoms by boundary 1009 mol->SetBoxDimension(&BoxLengths); 1010 *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " 1011 << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " 1012 << BoxLengths.x[2] << " with total volume of " << cellvolume << " " 1013 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl; 993 1014 } 994 1015 ; … … 1000 1021 Tesselation::Tesselation() 1001 1022 { 1002 1003 1004 1005 1023 PointsOnBoundaryCount = 0; 1024 LinesOnBoundaryCount = 0; 1025 TrianglesOnBoundaryCount = 0; 1026 TriangleFilesWritten = 0; 1006 1027 } 1007 1028 ; … … 1012 1033 Tesselation::~Tesselation() 1013 1034 { 1014 1015 1016 1017 1018 1019 1020 1021 1035 cout << Verbose(1) << "Free'ing TesselStruct ... " << endl; 1036 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { 1037 if (runner->second != NULL) { 1038 delete (runner->second); 1039 runner->second = NULL; 1040 } else 1041 cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl; 1042 } 1022 1043 } 1023 1044 ; … … 1031 1052 Tesselation::GuessStartingTriangle(ofstream *out) 1032 1053 { 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 //// listing distances1064 //*out << Verbose(1) << "Listing DistanceMMap:";1065 //for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {1066 //*out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";1067 //}1068 //*out << endl;1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1054 // 4b. create a starting triangle 1055 // 4b1. create all distances 1056 DistanceMultiMap DistanceMMap; 1057 double distance, tmp; 1058 Vector PlaneVector, TrialVector; 1059 PointMap::iterator A, B, C; // three nodes of the first triangle 1060 A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily 1061 1062 // with A chosen, take each pair B,C and sort 1063 if (A != PointsOnBoundary.end()) 1064 { 1065 B = A; 1066 B++; 1067 for (; B != PointsOnBoundary.end(); B++) 1068 { 1069 C = B; 1070 C++; 1071 for (; C != PointsOnBoundary.end(); C++) 1072 { 1073 tmp = A->second->node->x.DistanceSquared(&B->second->node->x); 1074 distance = tmp * tmp; 1075 tmp = A->second->node->x.DistanceSquared(&C->second->node->x); 1076 distance += tmp * tmp; 1077 tmp = B->second->node->x.DistanceSquared(&C->second->node->x); 1078 distance += tmp * tmp; 1079 DistanceMMap.insert(DistanceMultiMapPair(distance, pair< 1080 PointMap::iterator, PointMap::iterator> (B, C))); 1081 } 1082 } 1083 } 1084 // // listing distances 1085 // *out << Verbose(1) << "Listing DistanceMMap:"; 1086 // for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) { 1087 // *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")"; 1088 // } 1089 // *out << endl; 1090 // 4b2. pick three baselines forming a triangle 1091 // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate 1092 DistanceMultiMap::iterator baseline = DistanceMMap.begin(); 1093 for (; baseline != DistanceMMap.end(); baseline++) 1094 { 1095 // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate 1096 // 2. next, we have to check whether all points reside on only one side of the triangle 1097 // 3. construct plane vector 1098 PlaneVector.MakeNormalVector(&A->second->node->x, 1099 &baseline->second.first->second->node->x, 1100 &baseline->second.second->second->node->x); 1101 *out << Verbose(2) << "Plane vector of candidate triangle is "; 1102 PlaneVector.Output(out); 1103 *out << endl; 1104 // 4. loop over all points 1105 double sign = 0.; 1106 PointMap::iterator checker = PointsOnBoundary.begin(); 1107 for (; checker != PointsOnBoundary.end(); checker++) 1108 { 1109 // (neglecting A,B,C) 1110 if ((checker == A) || (checker == baseline->second.first) || (checker 1111 == baseline->second.second)) 1112 continue; 1113 // 4a. project onto plane vector 1114 TrialVector.CopyVector(&checker->second->node->x); 1115 TrialVector.SubtractVector(&A->second->node->x); 1116 distance = TrialVector.Projection(&PlaneVector); 1117 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok 1118 continue; 1119 *out << Verbose(3) << "Projection of " << checker->second->node->Name 1120 << " yields distance of " << distance << "." << endl; 1121 tmp = distance / fabs(distance); 1122 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle) 1123 if ((sign != 0) && (tmp != sign)) 1124 { 1125 // 4c. If so, break 4. loop and continue with next candidate in 1. loop 1126 *out << Verbose(2) << "Current candidates: " 1127 << A->second->node->Name << "," 1128 << baseline->second.first->second->node->Name << "," 1129 << baseline->second.second->second->node->Name << " leave " 1130 << checker->second->node->Name << " outside the convex hull." 1131 << endl; 1132 break; 1133 } 1134 else 1135 { // note the sign for later 1136 *out << Verbose(2) << "Current candidates: " 1137 << A->second->node->Name << "," 1138 << baseline->second.first->second->node->Name << "," 1139 << baseline->second.second->second->node->Name << " leave " 1140 << checker->second->node->Name << " inside the convex hull." 1141 << endl; 1142 sign = tmp; 1143 } 1144 // 4d. Check whether the point is inside the triangle (check distance to each node 1145 tmp = checker->second->node->x.DistanceSquared(&A->second->node->x); 1146 int innerpoint = 0; 1147 if ((tmp < A->second->node->x.DistanceSquared( 1148 &baseline->second.first->second->node->x)) && (tmp 1149 < A->second->node->x.DistanceSquared( 1150 &baseline->second.second->second->node->x))) 1151 innerpoint++; 1152 tmp = checker->second->node->x.DistanceSquared( 1153 &baseline->second.first->second->node->x); 1154 if ((tmp < baseline->second.first->second->node->x.DistanceSquared( 1155 &A->second->node->x)) && (tmp 1156 < baseline->second.first->second->node->x.DistanceSquared( 1157 &baseline->second.second->second->node->x))) 1158 innerpoint++; 1159 tmp = checker->second->node->x.DistanceSquared( 1160 &baseline->second.second->second->node->x); 1161 if ((tmp < baseline->second.second->second->node->x.DistanceSquared( 1162 &baseline->second.first->second->node->x)) && (tmp 1163 < baseline->second.second->second->node->x.DistanceSquared( 1164 &A->second->node->x))) 1165 innerpoint++; 1166 // 4e. If so, break 4. loop and continue with next candidate in 1. loop 1167 if (innerpoint == 3) 1168 break; 1169 } 1170 // 5. come this far, all on same side? Then break 1. loop and construct triangle 1171 if (checker == PointsOnBoundary.end()) 1172 { 1173 *out << "Looks like we have a candidate!" << endl; 1174 break; 1175 } 1176 } 1177 if (baseline != DistanceMMap.end()) 1178 { 1179 BPS[0] = baseline->second.first->second; 1180 BPS[1] = baseline->second.second->second; 1181 BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1182 BPS[0] = A->second; 1183 BPS[1] = baseline->second.second->second; 1184 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1185 BPS[0] = baseline->second.first->second; 1186 BPS[1] = A->second; 1187 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1188 1189 // 4b3. insert created triangle 1190 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1191 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1192 TrianglesOnBoundaryCount++; 1193 for (int i = 0; i < NDIM; i++) 1194 { 1195 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i])); 1196 LinesOnBoundaryCount++; 1197 } 1198 1199 *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl; 1200 } 1201 else 1202 { 1203 *out << Verbose(1) << "No starting triangle found." << endl; 1204 exit(255); 1205 } 1185 1206 } 1186 1207 ; … … 1191 1212 * -# if the lines contains to only one triangle 1192 1213 * -# We search all points in the boundary 1193 * 1194 * 1195 * 1196 * 1214 * -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors 1215 * -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to 1216 * baseline in triangle plane pointing out of the triangle and normal vector of new triangle) 1217 * -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount) 1197 1218 * \param *out output stream for debugging 1198 1219 * \param *configuration for IsAngstroem … … 1201 1222 void 1202 1223 Tesselation::TesselateOnBoundary(ofstream *out, config *configuration, 1203 1204 { 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 //if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())1292 //*out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;1293 //else1294 //*out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;1295 //if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())1296 //*out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;1297 //else1298 //*out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1224 molecule *mol) 1225 { 1226 bool flag; 1227 PointMap::iterator winner; 1228 class BoundaryPointSet *peak = NULL; 1229 double SmallestAngle, TempAngle; 1230 Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, 1231 PropagationVector; 1232 LineMap::iterator LineChecker[2]; 1233 do 1234 { 1235 flag = false; 1236 for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline 1237 != LinesOnBoundary.end(); baseline++) 1238 if (baseline->second->TrianglesCount == 1) 1239 { 1240 *out << Verbose(2) << "Current baseline is between " 1241 << *(baseline->second) << "." << endl; 1242 // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles) 1243 SmallestAngle = M_PI; 1244 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far 1245 // get peak point with respect to this base line's only triangle 1246 for (int i = 0; i < 3; i++) 1247 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) 1248 && (BTS->endpoints[i] != baseline->second->endpoints[1])) 1249 peak = BTS->endpoints[i]; 1250 *out << Verbose(3) << " and has peak " << *peak << "." << endl; 1251 // normal vector of triangle 1252 BTS->GetNormalVector(NormalVector); 1253 *out << Verbose(4) << "NormalVector of base triangle is "; 1254 NormalVector.Output(out); 1255 *out << endl; 1256 // offset to center of triangle 1257 CenterVector.Zero(); 1258 for (int i = 0; i < 3; i++) 1259 CenterVector.AddVector(&BTS->endpoints[i]->node->x); 1260 CenterVector.Scale(1. / 3.); 1261 *out << Verbose(4) << "CenterVector of base triangle is "; 1262 CenterVector.Output(out); 1263 *out << endl; 1264 // vector in propagation direction (out of triangle) 1265 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection) 1266 TempVector.CopyVector(&baseline->second->endpoints[0]->node->x); 1267 TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x); 1268 PropagationVector.MakeNormalVector(&TempVector, &NormalVector); 1269 TempVector.CopyVector(&CenterVector); 1270 TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center! 1271 //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl; 1272 if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline 1273 PropagationVector.Scale(-1.); 1274 *out << Verbose(4) << "PropagationVector of base triangle is "; 1275 PropagationVector.Output(out); 1276 *out << endl; 1277 winner = PointsOnBoundary.end(); 1278 for (PointMap::iterator target = PointsOnBoundary.begin(); target 1279 != PointsOnBoundary.end(); target++) 1280 if ((target->second != baseline->second->endpoints[0]) 1281 && (target->second != baseline->second->endpoints[1])) 1282 { // don't take the same endpoints 1283 *out << Verbose(3) << "Target point is " << *(target->second) 1284 << ":"; 1285 bool continueflag = true; 1286 1287 VirtualNormalVector.CopyVector( 1288 &baseline->second->endpoints[0]->node->x); 1289 VirtualNormalVector.AddVector( 1290 &baseline->second->endpoints[0]->node->x); 1291 VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line 1292 VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target 1293 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1294 continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees) 1295 if (!continueflag) 1296 { 1297 *out << Verbose(4) 1298 << "Angle between propagation direction and base line to " 1299 << *(target->second) << " is " << TempAngle 1300 << ", bad direction!" << endl; 1301 continue; 1302 } 1303 else 1304 *out << Verbose(4) 1305 << "Angle between propagation direction and base line to " 1306 << *(target->second) << " is " << TempAngle 1307 << ", good direction!" << endl; 1308 LineChecker[0] = baseline->second->endpoints[0]->lines.find( 1309 target->first); 1310 LineChecker[1] = baseline->second->endpoints[1]->lines.find( 1311 target->first); 1312 // if (LineChecker[0] != baseline->second->endpoints[0]->lines.end()) 1313 // *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl; 1314 // else 1315 // *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl; 1316 // if (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) 1317 // *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl; 1318 // else 1319 // *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl; 1320 // check first endpoint (if any connecting line goes to target or at least not more than 1) 1321 continueflag = continueflag && (((LineChecker[0] 1322 == baseline->second->endpoints[0]->lines.end()) 1323 || (LineChecker[0]->second->TrianglesCount == 1))); 1324 if (!continueflag) 1325 { 1326 *out << Verbose(4) << *(baseline->second->endpoints[0]) 1327 << " has line " << *(LineChecker[0]->second) 1328 << " to " << *(target->second) 1329 << " as endpoint with " 1330 << LineChecker[0]->second->TrianglesCount 1331 << " triangles." << endl; 1332 continue; 1333 } 1334 // check second endpoint (if any connecting line goes to target or at least not more than 1) 1335 continueflag = continueflag && (((LineChecker[1] 1336 == baseline->second->endpoints[1]->lines.end()) 1337 || (LineChecker[1]->second->TrianglesCount == 1))); 1338 if (!continueflag) 1339 { 1340 *out << Verbose(4) << *(baseline->second->endpoints[1]) 1341 << " has line " << *(LineChecker[1]->second) 1342 << " to " << *(target->second) 1343 << " as endpoint with " 1344 << LineChecker[1]->second->TrianglesCount 1345 << " triangles." << endl; 1346 continue; 1347 } 1348 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint) 1349 continueflag = continueflag && (!(((LineChecker[0] 1350 != baseline->second->endpoints[0]->lines.end()) 1351 && (LineChecker[1] 1352 != baseline->second->endpoints[1]->lines.end()) 1353 && (GetCommonEndpoint(LineChecker[0]->second, 1354 LineChecker[1]->second) == peak)))); 1355 if (!continueflag) 1356 { 1357 *out << Verbose(4) << "Current target is peak!" << endl; 1358 continue; 1359 } 1360 // in case NOT both were found 1361 if (continueflag) 1362 { // create virtually this triangle, get its normal vector, calculate angle 1363 flag = true; 1364 VirtualNormalVector.MakeNormalVector( 1365 &baseline->second->endpoints[0]->node->x, 1366 &baseline->second->endpoints[1]->node->x, 1367 &target->second->node->x); 1368 // make it always point inward 1369 if (baseline->second->endpoints[0]->node->x.Projection( 1370 &VirtualNormalVector) > 0) 1371 VirtualNormalVector.Scale(-1.); 1372 // calculate angle 1373 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1374 *out << Verbose(4) << "NormalVector is "; 1375 VirtualNormalVector.Output(out); 1376 *out << " and the angle is " << TempAngle << "." << endl; 1377 if (SmallestAngle > TempAngle) 1378 { // set to new possible winner 1379 SmallestAngle = TempAngle; 1380 winner = target; 1381 } 1382 } 1383 } 1384 // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle 1385 if (winner != PointsOnBoundary.end()) 1386 { 1387 *out << Verbose(2) << "Winning target point is " 1388 << *(winner->second) << " with angle " << SmallestAngle 1389 << "." << endl; 1390 // create the lins of not yet present 1391 BLS[0] = baseline->second; 1392 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles) 1393 LineChecker[0] = baseline->second->endpoints[0]->lines.find( 1394 winner->first); 1395 LineChecker[1] = baseline->second->endpoints[1]->lines.find( 1396 winner->first); 1397 if (LineChecker[0] 1398 == baseline->second->endpoints[0]->lines.end()) 1399 { // create 1400 BPS[0] = baseline->second->endpoints[0]; 1401 BPS[1] = winner->second; 1402 BLS[1] = new class BoundaryLineSet(BPS, 1403 LinesOnBoundaryCount); 1404 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, 1405 BLS[1])); 1406 LinesOnBoundaryCount++; 1407 } 1408 else 1409 BLS[1] = LineChecker[0]->second; 1410 if (LineChecker[1] 1411 == baseline->second->endpoints[1]->lines.end()) 1412 { // create 1413 BPS[0] = baseline->second->endpoints[1]; 1414 BPS[1] = winner->second; 1415 BLS[2] = new class BoundaryLineSet(BPS, 1416 LinesOnBoundaryCount); 1417 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, 1418 BLS[2])); 1419 LinesOnBoundaryCount++; 1420 } 1421 else 1422 BLS[2] = LineChecker[1]->second; 1423 BTS = new class BoundaryTriangleSet(BLS, 1424 TrianglesOnBoundaryCount); 1425 TrianglesOnBoundary.insert(TrianglePair( 1426 TrianglesOnBoundaryCount, BTS)); 1427 TrianglesOnBoundaryCount++; 1428 } 1429 else 1430 { 1431 *out << Verbose(1) 1432 << "I could not determine a winner for this baseline " 1433 << *(baseline->second) << "." << endl; 1434 } 1435 1436 // 5d. If the set of lines is not yet empty, go to 5. and continue 1437 } 1438 else 1439 *out << Verbose(2) << "Baseline candidate " << *(baseline->second) 1440 << " has a triangle count of " 1441 << baseline->second->TrianglesCount << "." << endl; 1442 } 1443 while (flag); 1423 1444 1424 1445 } … … 1431 1452 Tesselation::AddPoint(atom *Walker) 1432 1453 { 1433 1434 1435 1436 1437 1454 PointTestPair InsertUnique; 1455 BPS[0] = new class BoundaryPointSet(Walker); 1456 InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0])); 1457 if (InsertUnique.second) // if new point was not present before, increase counter 1458 PointsOnBoundaryCount++; 1438 1459 } 1439 1460 ; … … 1447 1468 Tesselation::AddTrianglePoint(atom* Candidate, int n) 1448 1469 { 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1470 PointTestPair InsertUnique; 1471 TPS[n] = new class BoundaryPointSet(Candidate); 1472 InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n])); 1473 if (InsertUnique.second) // if new point was not present before, increase counter 1474 { 1475 PointsOnBoundaryCount++; 1476 } 1477 else 1478 { 1479 delete TPS[n]; 1480 cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node) 1481 << " gibt's schon in der PointMap." << endl; 1482 TPS[n] = (InsertUnique.first)->second; 1483 } 1463 1484 } 1464 1485 ; 1465 1486 1466 1487 /** Function tries to add line from current Points in BPS to BoundaryLineSet. 1467 * If succes ful it raises the line count and inserts the new line into the BLS,1468 * if unsucces ful, it writes the line which had been present into the BLS, deleting the new constructed one.1488 * If successful it raises the line count and inserts the new line into the BLS, 1489 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one. 1469 1490 * @param *a first endpoint 1470 1491 * @param *b second endpoint 1471 1492 * @param n index of Tesselation::BLS giving the line with both endpoints 1472 1493 */ 1473 void 1474 Tesselation::AddTriangleLine(class BoundaryPointSet *a, 1475 class BoundaryPointSet *b, int n) 1476 { 1477 LineMap::iterator LineWalker; 1478 //cout << "Manually checking endpoints for line." << endl; 1479 if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr) 1480 //If a line is there, how do I recognize that beyond a shadow of a doubt? 1481 { 1482 //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl; 1483 1484 LineWalker = LinesOnBoundary.end(); 1485 LineWalker--; 1486 1487 while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr, 1488 b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max( 1489 a->node->nr, b->node->nr)) 1490 { 1491 //cout << Verbose(1) << "Looking for line which already exists"<< endl; 1492 LineWalker--; 1493 } 1494 BPS[0] = LineWalker->second->endpoints[0]; 1495 BPS[1] = LineWalker->second->endpoints[1]; 1496 BLS[n] = LineWalker->second; 1497 1498 } 1499 else 1500 { 1501 cout << Verbose(2) 1502 << "Adding line which has not been used before between " 1503 << *(a->node) << " and " << *(b->node) << "." << endl; 1504 BPS[0] = a; 1505 BPS[1] = b; 1506 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1507 1508 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n])); 1509 LinesOnBoundaryCount++; 1510 1511 } 1494 void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) { 1495 bool insertNewLine = true; 1496 1497 if (a->lines.find(b->node->nr) != a->lines.end()) { 1498 LineMap::iterator FindLine; 1499 pair<LineMap::iterator,LineMap::iterator> FindPair; 1500 FindPair = a->lines.equal_range(b->node->nr); 1501 1502 for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) { 1503 // If there is a line with less than two attached triangles, we don't need a new line. 1504 if (FindLine->second->TrianglesCount < 2) { 1505 insertNewLine = false; 1506 cout << Verbose(2) 1507 << "Using existing line " << *FindLine->second << endl; 1508 1509 BPS[0] = FindLine->second->endpoints[0]; 1510 BPS[1] = FindLine->second->endpoints[1]; 1511 BLS[n] = FindLine->second; 1512 1513 break; 1514 } 1515 } 1516 } 1517 1518 if (insertNewLine) { 1519 AlwaysAddTriangleLine(a, b, n); 1520 } 1512 1521 } 1513 1522 ; 1523 1524 /** 1525 * Adds lines from each of the current points in the BPS to BoundaryLineSet. 1526 * Raises the line count and inserts the new line into the BLS. 1527 * 1528 * @param *a first endpoint 1529 * @param *b second endpoint 1530 * @param n index of Tesselation::BLS giving the line with both endpoints 1531 */ 1532 void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) 1533 { 1534 cout << Verbose(2) 1535 << "Adding line which has not been used before between " 1536 << *(a->node) << " and " << *(b->node) << "." << endl; 1537 BPS[0] = a; 1538 BPS[1] = b; 1539 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1540 1541 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n])); 1542 LinesOnBoundaryCount++; 1543 }; 1514 1544 1515 1545 /** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm). … … 1520 1550 { 1521 1551 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1552 cout << Verbose(1) << "Adding triangle to its lines" << endl; 1553 int i = 0; 1554 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1555 TrianglesOnBoundaryCount++; 1556 1557 /* 1558 * this is apparently done when constructing triangle 1559 1560 for (i=0; i<3; i++) 1561 { 1562 BLS[i]->AddTriangle(BTS); 1563 } 1564 */ 1535 1565 } 1536 1566 ; … … 1538 1568 1539 1569 double det_get(gsl_matrix *A, int inPlace) { 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1570 /* 1571 inPlace = 1 => A is replaced with the LU decomposed copy. 1572 inPlace = 0 => A is retained, and a copy is used for LU. 1573 */ 1574 1575 double det; 1576 int signum; 1577 gsl_permutation *p = gsl_permutation_alloc(A->size1); 1578 gsl_matrix *tmpA; 1579 1580 if (inPlace) 1581 tmpA = A; 1582 else { 1583 gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2); 1584 gsl_matrix_memcpy(tmpA , A); 1585 } 1586 1587 1588 gsl_linalg_LU_decomp(tmpA , p , &signum); 1589 det = gsl_linalg_LU_det(tmpA , signum); 1590 gsl_permutation_free(p); 1591 if (! inPlace) 1592 gsl_matrix_free(tmpA); 1593 1594 return det; 1565 1595 }; 1566 1596 1567 1597 void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS) 1568 1598 { 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 center->x[0] =0.5 * m12/ m11;1604 1605 center->x[2] =0.5 * m14/ m11;1606 1607 1608 1609 1610 1599 gsl_matrix *A = gsl_matrix_calloc(3,3); 1600 double m11, m12, m13, m14; 1601 1602 for(int i=0;i<3;i++) { 1603 gsl_matrix_set(A, i, 0, a.x[i]); 1604 gsl_matrix_set(A, i, 1, b.x[i]); 1605 gsl_matrix_set(A, i, 2, c.x[i]); 1606 } 1607 m11 = det_get(A, 1); 1608 1609 for(int i=0;i<3;i++) { 1610 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]); 1611 gsl_matrix_set(A, i, 1, b.x[i]); 1612 gsl_matrix_set(A, i, 2, c.x[i]); 1613 } 1614 m12 = det_get(A, 1); 1615 1616 for(int i=0;i<3;i++) { 1617 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]); 1618 gsl_matrix_set(A, i, 1, a.x[i]); 1619 gsl_matrix_set(A, i, 2, c.x[i]); 1620 } 1621 m13 = det_get(A, 1); 1622 1623 for(int i=0;i<3;i++) { 1624 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]); 1625 gsl_matrix_set(A, i, 1, a.x[i]); 1626 gsl_matrix_set(A, i, 2, b.x[i]); 1627 } 1628 m14 = det_get(A, 1); 1629 1630 if (fabs(m11) < MYEPSILON) 1631 cerr << "ERROR: three points are colinear." << endl; 1632 1633 center->x[0] = 0.5 * m12/ m11; 1634 center->x[1] = -0.5 * m13/ m11; 1635 center->x[2] = 0.5 * m14/ m11; 1636 1637 if (fabs(a.Distance(center) - RADIUS) > MYEPSILON) 1638 cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl; 1639 1640 gsl_matrix_free(A); 1611 1641 }; 1612 1642 … … 1619 1649 * @param b vector second point of triangle 1620 1650 * @param c vector third point of triangle 1651 * @param *Umkreismittelpunkt new cneter point of circumference 1621 1652 * @param Direction vector indicates up/down 1622 1653 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle … … 1630 1661 */ 1631 1662 void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection, 1632 1633 { 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1663 double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius) 1664 { 1665 Vector TempNormal, helper; 1666 double Restradius; 1667 Vector OtherCenter; 1668 cout << Verbose(3) << "Begin of Get_center_of_sphere.\n"; 1669 Center->Zero(); 1670 helper.CopyVector(&a); 1671 helper.Scale(sin(2.*alpha)); 1672 Center->AddVector(&helper); 1673 helper.CopyVector(&b); 1674 helper.Scale(sin(2.*beta)); 1675 Center->AddVector(&helper); 1676 helper.CopyVector(&c); 1677 helper.Scale(sin(2.*gamma)); 1678 Center->AddVector(&helper); 1679 //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ; 1680 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma))); 1681 NewUmkreismittelpunkt->CopyVector(Center); 1682 cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n"; 1683 // Here we calculated center of circumscribing circle, using barycentric coordinates 1684 cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n"; 1685 1686 TempNormal.CopyVector(&a); 1687 TempNormal.SubtractVector(&b); 1688 helper.CopyVector(&a); 1689 helper.SubtractVector(&c); 1690 TempNormal.VectorProduct(&helper); 1691 if (fabs(HalfplaneIndicator) < MYEPSILON) 1692 { 1693 if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0)) 1694 { 1695 TempNormal.Scale(-1); 1696 } 1697 } 1698 else 1699 { 1700 if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0) 1701 { 1702 TempNormal.Scale(-1); 1703 } 1704 } 1705 1706 TempNormal.Normalize(); 1707 Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius); 1708 cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n"; 1709 TempNormal.Scale(Restradius); 1710 cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n"; 1711 1712 Center->AddVector(&TempNormal); 1713 cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n"; 1714 get_sphere(&OtherCenter, a, b, c, RADIUS); 1715 cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n"; 1716 cout << Verbose(3) << "End of Get_center_of_sphere.\n"; 1686 1717 }; 1718 1719 1687 1720 1688 1721 /** This recursive function finds a third point, to form a triangle with two given ones. 1689 1722 * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \ 1690 * 1691 * 1692 * 1693 * 1694 * 1695 * 1723 * supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \ 1724 * upon which we operate. 1725 * If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \ 1726 * direction and angle into Storage. 1727 * We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \ 1728 * with all neighbours of the candidate. 1696 1729 * @param a first point 1697 1730 * @param b second point … … 1710 1743 */ 1711 1744 void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, 1712 int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, 1713 atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol) 1714 { 1715 cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n"; 1716 cout << Verbose(3) << "Candidate is "<< *Candidate << endl; 1717 cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl; 1718 cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl; 1719 cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl; 1720 cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl; 1721 /* OldNormal is normal vector on the old triangle 1722 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal. 1723 * Chord points from b to a!!! 1724 */ 1725 Vector dif_a; //Vector from a to candidate 1726 Vector dif_b; //Vector from b to candidate 1727 Vector AngleCheck; 1728 Vector TempNormal, Umkreismittelpunkt; 1729 Vector Mittelpunkt; 1730 1731 double CurrentEpsilon = 0.1; 1732 double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance; 1733 double BallAngle, AlternativeSign; 1734 atom *Walker; // variable atom point 1735 1736 Vector NewUmkreismittelpunkt; 1737 1738 if (a != Candidate and b != Candidate and c != Candidate) { 1739 cout << Verbose(3) << "We have a unique candidate!" << endl; 1740 dif_a.CopyVector(&(a->x)); 1741 dif_a.SubtractVector(&(Candidate->x)); 1742 dif_b.CopyVector(&(b->x)); 1743 dif_b.SubtractVector(&(Candidate->x)); 1744 AngleCheck.CopyVector(&(Candidate->x)); 1745 AngleCheck.SubtractVector(&(a->x)); 1746 AngleCheck.ProjectOntoPlane(Chord); 1747 1748 SideA = dif_b.Norm(); 1749 SideB = dif_a.Norm(); 1750 SideC = Chord->Norm(); 1751 //Chord->Scale(-1); 1752 1753 alpha = Chord->Angle(&dif_a); 1754 beta = M_PI - Chord->Angle(&dif_b); 1755 gamma = dif_a.Angle(&dif_b); 1756 1757 cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl; 1758 1759 if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) { 1760 cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n"; 1761 cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl; 1762 } 1763 1764 if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) { 1765 Umkreisradius = SideA / 2.0 / sin(alpha); 1766 //cout << Umkreisradius << endl; 1767 //cout << SideB / 2.0 / sin(beta) << endl; 1768 //cout << SideC / 2.0 / sin(gamma) << endl; 1769 1770 if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points. 1771 cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl; 1772 cout << Verbose(2) << "Candidate is "<< *Candidate << endl; 1773 sign = AngleCheck.ScalarProduct(direction1); 1774 if (fabs(sign)<MYEPSILON) { 1775 if (AngleCheck.ScalarProduct(OldNormal)<0) { 1776 sign =0; 1777 AlternativeSign=1; 1778 } else { 1779 sign =0; 1780 AlternativeSign=-1; 1781 } 1782 } else { 1783 sign /= fabs(sign); 1784 } 1785 if (sign >= 0) { 1786 cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl; 1787 Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius); 1788 Mittelpunkt.SubtractVector(&ReferencePoint); 1789 cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl; 1790 BallAngle = Mittelpunkt.Angle(OldNormal); 1791 cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl; 1792 1793 //cout << "direction1 is " << *direction1 << "." << endl; 1794 //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl; 1795 //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl; 1796 1797 NewUmkreismittelpunkt.SubtractVector(&ReferencePoint); 1798 1799 if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) { 1800 if (Storage[0]< -1.5) { // first Candidate at all 1801 if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) { 1802 cout << Verbose(2) << "First good candidate is " << *Candidate << " with "; 1803 Opt_Candidate = Candidate; 1804 Storage[0] = sign; 1805 Storage[1] = AlternativeSign; 1806 Storage[2] = BallAngle; 1807 cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl; 1808 } else 1809 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl; 1810 } else { 1811 if ( Storage[2] > BallAngle) { 1812 if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) { 1813 cout << Verbose(2) << "Next better candidate is " << *Candidate << " with "; 1814 Opt_Candidate = Candidate; 1815 Storage[0] = sign; 1816 Storage[1] = AlternativeSign; 1817 Storage[2] = BallAngle; 1818 cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl; 1819 } else 1820 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl; 1821 } else { 1822 if (DEBUG) { 1823 cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl; 1824 } 1825 } 1826 } 1827 } else { 1828 if (DEBUG) { 1829 cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl; 1830 } 1831 } 1832 } else { 1833 if (DEBUG) { 1834 cout << Verbose(3) << *Candidate << " is not in search direction." << endl; 1835 } 1836 } 1837 } else { 1838 if (DEBUG) { 1839 cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl; 1840 } 1841 } 1842 } else { 1843 if (DEBUG) { 1844 cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl; 1845 } 1846 } 1847 } else { 1848 if (DEBUG) { 1849 cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl; 1850 } 1851 } 1852 1853 if (RecursionLevel < 5) { // Seven is the recursion level threshold. 1854 for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond 1855 Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate); 1856 if (Walker == Parent) { // don't go back the same bond 1857 continue; 1858 } else { 1859 Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again 1860 } 1861 } 1862 } 1863 cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n"; 1745 int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, 1746 atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol) 1747 { 1748 cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n"; 1749 cout << Verbose(3) << "Candidate is "<< *Candidate << endl; 1750 cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl; 1751 cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl; 1752 cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl; 1753 cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl; 1754 /* OldNormal is normal vector on the old triangle 1755 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal. 1756 * Chord points from b to a!!! 1757 */ 1758 Vector dif_a; //Vector from a to candidate 1759 Vector dif_b; //Vector from b to candidate 1760 Vector AngleCheck; 1761 Vector TempNormal, Umkreismittelpunkt; 1762 Vector Mittelpunkt; 1763 1764 double CurrentEpsilon = 0.1; 1765 double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance; 1766 double BallAngle, AlternativeSign; 1767 atom *Walker; // variable atom point 1768 1769 Vector NewUmkreismittelpunkt; 1770 1771 if (a != Candidate and b != Candidate and c != Candidate) { 1772 cout << Verbose(3) << "We have a unique candidate!" << endl; 1773 dif_a.CopyVector(&(a->x)); 1774 dif_a.SubtractVector(&(Candidate->x)); 1775 dif_b.CopyVector(&(b->x)); 1776 dif_b.SubtractVector(&(Candidate->x)); 1777 AngleCheck.CopyVector(&(Candidate->x)); 1778 AngleCheck.SubtractVector(&(a->x)); 1779 AngleCheck.ProjectOntoPlane(Chord); 1780 1781 SideA = dif_b.Norm(); 1782 SideB = dif_a.Norm(); 1783 SideC = Chord->Norm(); 1784 //Chord->Scale(-1); 1785 1786 alpha = Chord->Angle(&dif_a); 1787 beta = M_PI - Chord->Angle(&dif_b); 1788 gamma = dif_a.Angle(&dif_b); 1789 1790 cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl; 1791 1792 if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) { 1793 cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n"; 1794 cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl; 1795 } 1796 1797 if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) { 1798 Umkreisradius = SideA / 2.0 / sin(alpha); 1799 //cout << Umkreisradius << endl; 1800 //cout << SideB / 2.0 / sin(beta) << endl; 1801 //cout << SideC / 2.0 / sin(gamma) << endl; 1802 1803 if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points. 1804 cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl; 1805 cout << Verbose(2) << "Candidate is "<< *Candidate << endl; 1806 sign = AngleCheck.ScalarProduct(direction1); 1807 if (fabs(sign)<MYEPSILON) { 1808 if (AngleCheck.ScalarProduct(OldNormal)<0) { 1809 sign =0; 1810 AlternativeSign=1; 1811 } else { 1812 sign =0; 1813 AlternativeSign=-1; 1814 } 1815 } else { 1816 sign /= fabs(sign); 1817 } 1818 if (sign >= 0) { 1819 cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl; 1820 1821 Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius); 1822 1823 Mittelpunkt.SubtractVector(&ReferencePoint); 1824 cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl; 1825 1826 BallAngle = Mittelpunkt.Angle(OldNormal); 1827 cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl; 1828 1829 //cout << "direction1 is " << *direction1 << "." << endl; 1830 //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl; 1831 1832 //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl; 1833 1834 NewUmkreismittelpunkt.SubtractVector(&ReferencePoint); 1835 1836 if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) { 1837 if (Storage[0]< -1.5) { // first Candidate at all 1838 if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) { 1839 cout << Verbose(2) << "First good candidate is " << *Candidate << " with "; 1840 Opt_Candidate = Candidate; 1841 Storage[0] = sign; 1842 Storage[1] = AlternativeSign; 1843 Storage[2] = BallAngle; 1844 cout << " angle " << Storage[2] << " and Up/Down " 1845 << Storage[0] << endl; 1846 } else 1847 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl; 1848 } else { 1849 if ( Storage[2] > BallAngle) { 1850 if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) { 1851 cout << Verbose(2) << "Next better candidate is " << *Candidate << " with "; 1852 Opt_Candidate = Candidate; 1853 Storage[0] = sign; 1854 Storage[1] = AlternativeSign; 1855 Storage[2] = BallAngle; 1856 cout << " angle " << Storage[2] << " and Up/Down " 1857 << Storage[0] << endl; 1858 } else 1859 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl; 1860 } else { 1861 if (DEBUG) { 1862 cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl; 1863 } 1864 } 1865 } 1866 } else { 1867 if (DEBUG) { 1868 cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl; 1869 } 1870 } 1871 } else { 1872 if (DEBUG) { 1873 cout << Verbose(3) << *Candidate << " is not in search direction." << endl; 1874 } 1875 } 1876 } else { 1877 if (DEBUG) { 1878 cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl; 1879 } 1880 } 1881 } else { 1882 if (DEBUG) { 1883 cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl; 1884 } 1885 } 1886 } else { 1887 if (DEBUG) { 1888 cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl; 1889 } 1890 } 1891 1892 if (RecursionLevel < 5) { // Seven is the recursion level threshold. 1893 for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond 1894 Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate); 1895 if (Walker == Parent) { // don't go back the same bond 1896 continue; 1897 } else { 1898 Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again 1899 } 1900 } 1901 } 1902 cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n"; 1864 1903 } 1865 1904 ; … … 1874 1913 void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c) 1875 1914 { 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1915 Vector helper; 1916 double alpha, beta, gamma; 1917 Vector SideA, SideB, SideC; 1918 SideA.CopyVector(b); 1919 SideA.SubtractVector(c); 1920 SideB.CopyVector(c); 1921 SideB.SubtractVector(a); 1922 SideC.CopyVector(a); 1923 SideC.SubtractVector(b); 1924 alpha = M_PI - SideB.Angle(&SideC); 1925 beta = M_PI - SideC.Angle(&SideA); 1926 gamma = M_PI - SideA.Angle(&SideB); 1927 cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl; 1928 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) 1929 cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl; 1930 1931 Center->Zero(); 1932 helper.CopyVector(a); 1933 helper.Scale(sin(2.*alpha)); 1934 Center->AddVector(&helper); 1935 helper.CopyVector(b); 1936 helper.Scale(sin(2.*beta)); 1937 Center->AddVector(&helper); 1938 helper.CopyVector(c); 1939 helper.Scale(sin(2.*gamma)); 1940 Center->AddVector(&helper); 1941 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma))); 1903 1942 }; 1904 1943 … … 1918 1957 double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection) 1919 1958 { 1920 1921 1922 1923 1924 1925 1926 cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))<< "!" << endl;1927 1928 1929 1930 1931 1932 1933 1934 1935 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;1946 1947 1959 Vector helper; 1960 double radius, alpha; 1961 1962 helper.CopyVector(&NewSphereCenter); 1963 // test whether new center is on the parameter circle's plane 1964 if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 1965 cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 1966 helper.ProjectOntoPlane(&CirclePlaneNormal); 1967 } 1968 radius = helper.ScalarProduct(&helper); 1969 // test whether the new center vector has length of CircleRadius 1970 if (fabs(radius - CircleRadius) > HULLEPSILON) 1971 cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 1972 alpha = helper.Angle(&OldSphereCenter); 1973 // make the angle unique by checking the halfplanes/search direction 1974 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals 1975 alpha = 2.*M_PI - alpha; 1976 cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl; 1977 radius = helper.Distance(&OldSphereCenter); 1978 helper.ProjectOntoPlane(&NormalVector); 1979 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles 1980 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) { 1981 cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl; 1982 return alpha; 1983 } else { 1984 cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl; 1985 return 2.*M_PI; 1986 } 1948 1987 }; 1949 1988 … … 1980 2019 // void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC) 1981 2020 // { 1982 // 1983 // Vector CircleCenter;// center of the circle, i.e. of the band of sphere's centers1984 // 1985 // Vector OldSphereCenter;// center of the sphere defined by the three points of BaseTriangle1986 // Vector NewSphereCenter;// center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility1987 // Vector OtherNewSphereCenter;// center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility1988 // Vector NewNormalVector;// normal vector of the Candidate's triangle1989 // Vector SearchDirection;// vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)1990 // 1991 // 1992 // 1993 // 1994 // 1995 // 1996 // 1997 // 2021 // atom *Walker = NULL; 2022 // Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2023 // Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 2024 // Vector OldSphereCenter; // center of the sphere defined by the three points of BaseTriangle 2025 // Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 2026 // Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 2027 // Vector NewNormalVector; // normal vector of the Candidate's triangle 2028 // Vector SearchDirection; // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate) 2029 // Vector helper; 2030 // LinkedAtoms *List = NULL; 2031 // double CircleRadius; // radius of this circle 2032 // double radius; 2033 // double alpha, Otheralpha; // angles (i.e. parameter for the circle). 2034 // double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle 2035 // int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2036 // atom *Candidate = NULL; 1998 2037 // 1999 // 2038 // cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl; 2000 2039 // 2001 // 2040 // cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl; 2002 2041 // 2003 // 2004 // 2005 // 2006 // 2042 // // construct center of circle 2043 // CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x)); 2044 // CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x); 2045 // CircleCenter.Scale(0.5); 2007 2046 // 2008 // 2009 // 2010 // 2047 // // construct normal vector of circle 2048 // CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x); 2049 // CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x); 2011 2050 // 2012 // 2013 // 2014 // 2015 // 2016 // 2017 // 2051 // // calculate squared radius of circle 2052 // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2053 // if (radius/4. < RADIUS*RADIUS) { 2054 // CircleRadius = RADIUS*RADIUS - radius/4.; 2055 // CirclePlaneNormal.Normalize(); 2056 // cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2018 2057 // 2019 // 2020 // 2021 // helper.CopyVector(&BaseTriangle->NormalVector);// normal vector ensures that this is correct center of the two possible ones2022 // 2023 // 2024 // 2025 // 2026 // 2058 // // construct old center 2059 // GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x)); 2060 // helper.CopyVector(&BaseTriangle->NormalVector); // normal vector ensures that this is correct center of the two possible ones 2061 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter); 2062 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2063 // OldSphereCenter.AddVector(&helper); 2064 // OldSphereCenter.SubtractVector(&CircleCenter); 2065 // cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2027 2066 // 2028 // 2029 // 2030 // 2031 // 2032 // 2033 // 2034 // 2067 // // test whether old center is on the band's plane 2068 // if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 2069 // cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 2070 // OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 2071 // } 2072 // radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2073 // if (fabs(radius - CircleRadius) < HULLEPSILON) { 2035 2074 // 2036 // 2037 // 2038 // 2039 // for(int i=0;i<3;i++)// just take next different endpoint2040 // 2041 // 2042 // 2043 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!2044 // 2045 // 2046 // 2047 // 2048 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {// rotated the wrong way!2049 // 2050 // 2075 // // construct SearchDirection 2076 // SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal); 2077 // helper.CopyVector(&BaseLine->endpoints[0]->node->x); 2078 // for(int i=0;i<3;i++) // just take next different endpoint 2079 // if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) { 2080 // helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x); 2081 // } 2082 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards! 2083 // SearchDirection.Scale(-1.); 2084 // SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2085 // SearchDirection.Normalize(); 2086 // cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2087 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2088 // cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2089 // } 2051 2090 // 2052 // if (LC->SetIndexToVector(&CircleCenter)) {// get cell for the starting atom2053 // 2054 // 2055 // 2056 // 2057 // 2058 // 2059 // 2060 // 2061 // 2062 // 2063 // 2064 // 2065 // 2066 // 2067 // 2068 // 2069 // 2070 // 2071 // 2072 // 2073 // 2074 // 2075 // 2091 // if (LC->SetIndexToVector(&CircleCenter)) { // get cell for the starting atom 2092 // for(int i=0;i<NDIM;i++) // store indices of this cell 2093 // N[i] = LC->n[i]; 2094 // cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2095 // } else { 2096 // cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; 2097 // return; 2098 // } 2099 // // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2100 // cout << Verbose(2) << "LC Intervals:"; 2101 // for (int i=0;i<NDIM;i++) { 2102 // Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2103 // Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2104 // cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2105 // } 2106 // cout << endl; 2107 // for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2108 // for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2109 // for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2110 // List = LC->GetCurrentCell(); 2111 // cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2112 // if (List != NULL) { 2113 // for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2114 // Candidate = (*Runner); 2076 2115 // 2077 // 2078 // 2079 // 2116 // // check for three unique points 2117 // if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) { 2118 // cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl; 2080 2119 // 2081 // 2082 // 2083 // 2120 // // construct both new centers 2121 // GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)); 2122 // OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2084 2123 // 2085 // 2086 // 2087 // 2088 // 2089 // 2090 // 2091 // 2092 // 2093 // 2094 // 2124 // if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) { 2125 // helper.CopyVector(&NewNormalVector); 2126 // cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2127 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter); 2128 // if (radius < RADIUS*RADIUS) { 2129 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2130 // cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl; 2131 // NewSphereCenter.AddVector(&helper); 2132 // NewSphereCenter.SubtractVector(&CircleCenter); 2133 // cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2095 2134 // 2096 // 2097 // 2098 // 2099 // 2135 // helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction 2136 // OtherNewSphereCenter.AddVector(&helper); 2137 // OtherNewSphereCenter.SubtractVector(&CircleCenter); 2138 // cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2100 2139 // 2101 // 2102 // 2103 // 2104 // 2105 // 2106 // 2107 // 2108 // 2109 // 2110 // 2111 // 2112 // 2113 // 2114 // 2115 // 2116 // 2117 // 2118 // 2140 // // check both possible centers 2141 // alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection); 2142 // Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection); 2143 // alpha = min(alpha, Otheralpha); 2144 // if (*ShortestAngle > alpha) { 2145 // OptCandidate = Candidate; 2146 // *ShortestAngle = alpha; 2147 // if (alpha != Otheralpha) 2148 // OptCandidateCenter->CopyVector(&NewSphereCenter); 2149 // else 2150 // OptCandidateCenter->CopyVector(&OtherNewSphereCenter); 2151 // cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl; 2152 // } else { 2153 // if (OptCandidate != NULL) 2154 // cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2155 // else 2156 // cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2157 // } 2119 2158 // 2120 // 2121 // 2122 // 2123 // 2124 // 2125 // 2126 // 2127 // 2128 // 2129 // 2130 // 2131 // 2132 // 2133 // 2134 // 2135 // 2136 // 2137 // 2159 // } else { 2160 // cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl; 2161 // } 2162 // } else { 2163 // cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2164 // } 2165 // } else { 2166 // cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl; 2167 // } 2168 // } 2169 // } 2170 // } 2171 // } else { 2172 // cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 2173 // } 2174 // } else { 2175 // cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl; 2176 // } 2138 2177 // 2139 // 2178 // cout << Verbose(1) << "End of Find_next_suitable_point" << endl; 2140 2179 // }; 2141 2180 … … 2147 2186 * \param *out output stream for debugging 2148 2187 * \param *Candidates endpoints of the triangle candidate 2149 * \return false - triangle invalid due to edge criteria, true - triangle may be added. 2188 * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two 2189 * triangles exist which is the maximum for three points 2150 2190 */ 2151 bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) { 2152 LineMap::iterator FindLine; 2153 PointMap::iterator FindPoint; 2154 bool Present[3]; 2155 2156 *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl; 2157 for (int i=0;i<3;i++) { // check through all endpoints 2158 FindPoint = PointsOnBoundary.find(Candidates[i]->nr); 2159 if (FindPoint != PointsOnBoundary.end()) 2160 TPS[i] = FindPoint->second; 2161 else 2162 TPS[i] = NULL; 2163 } 2164 2165 // check lines 2166 for (int i=0;i<3;i++) 2167 if (TPS[i] != NULL) 2168 for (int j=i;j<3;j++) 2169 if (TPS[j] != NULL) { 2170 FindLine = TPS[i]->lines.find(TPS[j]->node->nr); 2171 if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) { 2172 *out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->second->TrianglesCount << " triangles attached." << endl; 2173 *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 2174 return false; 2175 } 2176 } 2177 *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 2178 return true; 2191 int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) { 2192 LineMap::iterator FindLine; 2193 PointMap::iterator FindPoint; 2194 TriangleMap::iterator FindTriangle; 2195 int adjacentTriangleCount = 0; 2196 class BoundaryPointSet *Points[3]; 2197 2198 *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl; 2199 // builds a triangle point set (Points) of the end points 2200 for (int i = 0; i < 3; i++) { 2201 FindPoint = PointsOnBoundary.find(Candidates[i]->nr); 2202 if (FindPoint != PointsOnBoundary.end()) { 2203 Points[i] = FindPoint->second; 2204 } else { 2205 Points[i] = NULL; 2206 } 2207 } 2208 2209 // checks lines between the points in the Points for their adjacent triangles 2210 for (int i = 0; i < 3; i++) { 2211 if (Points[i] != NULL) { 2212 for (int j = i; j < 3; j++) { 2213 if (Points[j] != NULL) { 2214 FindLine = Points[i]->lines.find(Points[j]->node->nr); 2215 if (FindLine != Points[i]->lines.end()) { 2216 for (; FindLine->first == Points[j]->node->nr; FindLine++) { 2217 FindTriangle = FindLine->second->triangles.begin(); 2218 for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 2219 if (( 2220 (FindTriangle->second->endpoints[0] == Points[0]) 2221 || (FindTriangle->second->endpoints[0] == Points[1]) 2222 || (FindTriangle->second->endpoints[0] == Points[2]) 2223 ) && ( 2224 (FindTriangle->second->endpoints[1] == Points[0]) 2225 || (FindTriangle->second->endpoints[1] == Points[1]) 2226 || (FindTriangle->second->endpoints[1] == Points[2]) 2227 ) && ( 2228 (FindTriangle->second->endpoints[2] == Points[0]) 2229 || (FindTriangle->second->endpoints[2] == Points[1]) 2230 || (FindTriangle->second->endpoints[2] == Points[2]) 2231 ) 2232 ) { 2233 adjacentTriangleCount++; 2234 } 2235 } 2236 } 2237 // Only one of the triangle lines must be considered for the triangle count. 2238 *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 2239 return adjacentTriangleCount; 2240 2241 } 2242 } 2243 } 2244 } 2245 } 2246 2247 *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 2248 return adjacentTriangleCount; 2179 2249 }; 2180 2250 … … 2206 2276 * @param BaseLine BoundaryLineSet with the current base line 2207 2277 * @param ThirdNode third atom to avoid in search 2208 * @param OptCandidate candidate reference on return 2209 * @param OptCandidateCenter candidate's sphere center on return 2278 * @param candidates list of equally good candidates to return 2210 2279 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate 2211 2280 * @param RADIUS radius of sphere 2212 2281 * @param *LC LinkedCell structure with neighbouring atoms 2213 2282 */ 2214 void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC) 2215 { 2216 atom *Walker = NULL; 2217 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2218 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 2219 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 2220 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 2221 Vector NewNormalVector; // normal vector of the Candidate's triangle 2222 Vector helper; 2223 LinkedAtoms *List = NULL; 2224 double CircleRadius; // radius of this circle 2225 double radius; 2226 double alpha, Otheralpha; // angles (i.e. parameter for the circle). 2227 double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle 2228 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2229 atom *Candidate = NULL; 2230 2231 cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl; 2232 2233 cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2234 2235 // construct center of circle 2236 CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x)); 2237 CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x); 2238 CircleCenter.Scale(0.5); 2239 2240 // construct normal vector of circle 2241 CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x); 2242 CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x); 2243 2244 // calculate squared radius of circle 2245 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2246 if (radius/4. < RADIUS*RADIUS) { 2247 CircleRadius = RADIUS*RADIUS - radius/4.; 2248 CirclePlaneNormal.Normalize(); 2249 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2250 2251 // test whether old center is on the band's plane 2252 if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 2253 cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 2254 OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 2255 } 2256 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2257 if (fabs(radius - CircleRadius) < HULLEPSILON) { 2258 2259 // check SearchDirection 2260 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2261 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2262 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; 2263 } 2264 // get cell for the starting atom 2265 if (LC->SetIndexToVector(&CircleCenter)) { 2266 for(int i=0;i<NDIM;i++) // store indices of this cell 2267 N[i] = LC->n[i]; 2268 cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2269 } else { 2270 cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; 2271 return; 2272 } 2273 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2274 cout << Verbose(2) << "LC Intervals:"; 2275 for (int i=0;i<NDIM;i++) { 2276 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2277 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2278 cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2279 } 2280 cout << endl; 2281 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2282 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2283 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2284 List = LC->GetCurrentCell(); 2285 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2286 if (List != NULL) { 2287 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2288 Candidate = (*Runner); 2289 2290 // check for three unique points 2291 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) { 2292 cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl; 2293 2294 // construct both new centers 2295 GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)); 2296 OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2297 2298 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) { 2299 helper.CopyVector(&NewNormalVector); 2300 cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2301 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter); 2302 if (radius < RADIUS*RADIUS) { 2303 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2304 cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl; 2305 NewSphereCenter.AddVector(&helper); 2306 NewSphereCenter.SubtractVector(&CircleCenter); 2307 cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2308 2309 helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction 2310 OtherNewSphereCenter.AddVector(&helper); 2311 OtherNewSphereCenter.SubtractVector(&CircleCenter); 2312 cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2313 2314 // check both possible centers 2315 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2316 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2317 alpha = min(alpha, Otheralpha); 2318 if (*ShortestAngle > alpha) { 2319 OptCandidate = Candidate; 2320 *ShortestAngle = alpha; 2321 if (alpha != Otheralpha) 2322 OptCandidateCenter->CopyVector(&NewSphereCenter); 2323 else 2324 OptCandidateCenter->CopyVector(&OtherNewSphereCenter); 2325 cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl; 2326 } else { 2327 if (OptCandidate != NULL) 2328 cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2329 else 2330 cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2331 } 2332 2333 } else { 2334 cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl; 2335 } 2336 } else { 2337 cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2338 } 2339 } else { 2340 if (ThirdNode != NULL) 2341 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2342 else 2343 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl; 2344 } 2345 } 2346 } 2347 } 2348 } else { 2349 cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 2350 } 2283 void Find_third_point_for_Tesselation( 2284 Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, 2285 class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates, 2286 double *ShortestAngle, const double RADIUS, LinkedCell *LC 2287 ) { 2288 atom *Walker = NULL; 2289 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2290 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 2291 Vector SphereCenter; 2292 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 2293 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 2294 Vector NewNormalVector; // normal vector of the Candidate's triangle 2295 Vector helper, OptCandidateCenter, OtherOptCandidateCenter; 2296 LinkedAtoms *List = NULL; 2297 double CircleRadius; // radius of this circle 2298 double radius; 2299 double alpha, Otheralpha; // angles (i.e. parameter for the circle). 2300 double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle 2301 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2302 atom *Candidate = NULL; 2303 CandidateForTesselation *optCandidate; 2304 2305 cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl; 2306 2307 cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2308 2309 // construct center of circle 2310 CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x)); 2311 CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x); 2312 CircleCenter.Scale(0.5); 2313 2314 // construct normal vector of circle 2315 CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x); 2316 CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x); 2317 2318 // calculate squared radius atom *ThirdNode,f circle 2319 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2320 if (radius/4. < RADIUS*RADIUS) { 2321 CircleRadius = RADIUS*RADIUS - radius/4.; 2322 CirclePlaneNormal.Normalize(); 2323 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2324 2325 // test whether old center is on the band's plane 2326 if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 2327 cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 2328 OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 2329 } 2330 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2331 if (fabs(radius - CircleRadius) < HULLEPSILON) { 2332 2333 // check SearchDirection 2334 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2335 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2336 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; 2337 } 2338 2339 // get cell for the starting atom 2340 if (LC->SetIndexToVector(&CircleCenter)) { 2341 for(int i=0;i<NDIM;i++) // store indices of this cell 2342 N[i] = LC->n[i]; 2343 cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2344 } else { 2345 cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; 2346 return; 2347 } 2348 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2349 cout << Verbose(2) << "LC Intervals:"; 2350 for (int i=0;i<NDIM;i++) { 2351 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2352 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2353 cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2354 } 2355 cout << endl; 2356 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2357 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2358 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2359 List = LC->GetCurrentCell(); 2360 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2361 if (List != NULL) { 2362 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2363 Candidate = (*Runner); 2364 2365 // check for three unique points 2366 cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl; 2367 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){ 2368 2369 // construct both new centers 2370 GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)); 2371 OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2372 2373 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) 2374 && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON) 2375 ) { 2376 helper.CopyVector(&NewNormalVector); 2377 cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2378 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter); 2379 if (radius < RADIUS*RADIUS) { 2380 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2381 cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl; 2382 NewSphereCenter.AddVector(&helper); 2383 NewSphereCenter.SubtractVector(&CircleCenter); 2384 cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2385 2386 // OtherNewSphereCenter is created by the same vector just in the other direction 2387 helper.Scale(-1.); 2388 OtherNewSphereCenter.AddVector(&helper); 2389 OtherNewSphereCenter.SubtractVector(&CircleCenter); 2390 cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2391 2392 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2393 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2394 alpha = min(alpha, Otheralpha); 2395 // if there is a better candidate, drop the current list and add the new candidate 2396 // otherwise ignore the new candidate and keep the list 2397 if (*ShortestAngle > (alpha - HULLEPSILON)) { 2398 optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter); 2399 if (fabs(alpha - Otheralpha) > MYEPSILON) { 2400 optCandidate->OptCenter.CopyVector(&NewSphereCenter); 2401 optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter); 2402 } else { 2403 optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter); 2404 optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter); 2405 } 2406 // if there is an equal candidate, add it to the list without clearing the list 2407 if ((*ShortestAngle - HULLEPSILON) < alpha) { 2408 candidates->push_back(optCandidate); 2409 cout << Verbose(1) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with " 2410 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2411 } else { 2412 candidates->clear(); 2413 candidates->push_back(optCandidate); 2414 cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with " 2415 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2416 } 2417 *ShortestAngle = alpha; 2418 cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl; 2419 } else { 2420 if ((optCandidate != NULL) && (optCandidate->point != NULL)) 2421 cout << Verbose(1) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2422 else 2423 cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2424 } 2425 2426 } else { 2427 cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl; 2428 } 2429 } else { 2430 cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2431 } 2432 } else { 2433 if (ThirdNode != NULL) 2434 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2435 else 2436 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl; 2437 } 2438 } 2439 } 2440 } 2441 } else { 2442 cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 2443 } 2444 } else { 2445 if (ThirdNode != NULL) 2446 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 2447 else 2448 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2449 } 2450 2451 cout << Verbose(1) << "INFO: Sorting candidate list ..." << endl; 2452 if (candidates->size() > 1) { 2453 candidates->unique(); 2454 candidates->sort(sortCandidates); 2455 } 2456 2457 cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl; 2458 }; 2459 2460 /** 2461 * Finds the preferable out of two third-point candidates with equal angles. 2462 * 2463 * @param Candidate - this and the second parameter are evaluated 2464 * @param OptCandidate - this and the second parameter are evaluated 2465 * @param current base line 2466 * @param third node of the base triangle 2467 * @param tesselation object 2468 * 2469 * @return true if Candidate should be taken, false if OptCandidate should be kept 2470 */ 2471 bool Choose_preferable_third_point( 2472 atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, 2473 atom *ThirdNode, Tesselation *Tess 2474 ) { 2475 bool takeNewCandidate; 2476 2477 ofstream *out = new ofstream(); 2478 atom *Atoms[3]; 2479 bool optCandidateAndBaseLineFormTriangle = (ThirdNode != NULL) && (OptCandidate == ThirdNode); 2480 bool candidateAndBaseLineFormTriangle = (ThirdNode != NULL) && (Candidate == ThirdNode); 2481 Atoms[0] = Candidate; 2482 Atoms[1] = OptCandidate; 2483 Atoms[2] = BaseLine->endpoints[0]->node; 2484 bool candidatesAndBaseLineNode0FormTriangle = (Tess->CheckPresenceOfTriangle(out, Atoms) > 0); 2485 Atoms[0] = Candidate; 2486 Atoms[1] = OptCandidate; 2487 Atoms[2] = BaseLine->endpoints[1]->node; 2488 bool candidatesAndBaseLineNode1FormTriangle = (Tess->CheckPresenceOfTriangle(out, Atoms) > 0); 2489 Vector halfBaseLine; 2490 halfBaseLine.CopyVector(&BaseLine->endpoints[0]->node->x); 2491 halfBaseLine.AddVector(&BaseLine->endpoints[1]->node->x); 2492 halfBaseLine.Scale(0.5); 2493 2494 if (optCandidateAndBaseLineFormTriangle) { 2495 takeNewCandidate = (!existsIntersection(Candidate->x, halfBaseLine, OptCandidate->x, BaseLine->endpoints[0]->node->x) 2496 && !existsIntersection(Candidate->x, halfBaseLine, OptCandidate->x, BaseLine->endpoints[1]->node->x)); 2497 } else if (candidateAndBaseLineFormTriangle) { 2498 takeNewCandidate = (existsIntersection(OptCandidate->x, halfBaseLine, Candidate->x, BaseLine->endpoints[0]->node->x) 2499 || existsIntersection(OptCandidate->x, halfBaseLine, Candidate->x, BaseLine->endpoints[1]->node->x)); 2500 } else if (candidatesAndBaseLineNode0FormTriangle) { 2501 takeNewCandidate = !existsIntersection(OptCandidate->x, BaseLine->endpoints[0]->node->x, Candidate->x, halfBaseLine); 2502 } else if (candidatesAndBaseLineNode1FormTriangle) { 2503 takeNewCandidate = !existsIntersection(OptCandidate->x, BaseLine->endpoints[1]->node->x, Candidate->x, halfBaseLine); 2351 2504 } else { 2352 if (ThirdNode != NULL) 2353 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 2354 else 2355 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2505 takeNewCandidate = (ThirdNode == NULL) 2506 || ((!existsIntersection(Candidate->x, halfBaseLine, ThirdNode->x, BaseLine->endpoints[0]->node->x) 2507 && !existsIntersection(Candidate->x, halfBaseLine, ThirdNode->x, BaseLine->endpoints[1]->node->x))); 2356 2508 } 2357 2509 2358 cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;2510 return takeNewCandidate; 2359 2511 }; 2512 2513 2514 struct Intersection { 2515 Vector x1; 2516 Vector x2; 2517 Vector x3; 2518 Vector x4; 2519 }; 2520 2521 /** 2522 * Intersection calculation function. 2523 * 2524 * @param x to find the result for 2525 * @param function parameter 2526 */ 2527 double MinIntersectDistance(const gsl_vector * x, void *params) { 2528 double retval = 0; 2529 struct Intersection *I = (struct Intersection *)params; 2530 Vector intersection; 2531 Vector SideA,SideB,HeightA, HeightB; 2532 for (int i=0;i<NDIM;i++) 2533 intersection.x[i] = gsl_vector_get(x, i); 2534 2535 SideA.CopyVector(&(I->x1)); 2536 SideA.SubtractVector(&I->x2); 2537 HeightA.CopyVector(&intersection); 2538 HeightA.SubtractVector(&I->x1); 2539 HeightA.ProjectOntoPlane(&SideA); 2540 2541 SideB.CopyVector(&I->x3); 2542 SideB.SubtractVector(&I->x4); 2543 HeightB.CopyVector(&intersection); 2544 HeightB.SubtractVector(&I->x3); 2545 HeightB.ProjectOntoPlane(&SideB); 2546 2547 retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB); 2548 //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl; 2549 2550 return retval; 2551 }; 2552 2553 2554 /** 2555 * Calculates whether there is an intersection between two lines. The first line 2556 * always goes through point 1 and point 2 and the second line is given by the 2557 * connection between point 4 and point 5. 2558 * 2559 * @param point 1 of line 1 2560 * @param point 2 of line 1 2561 * @param point 1 of line 2 2562 * @param point 2 of line 2 2563 * 2564 * @return true if there is an intersection between the given lines, false otherwise 2565 */ 2566 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) { 2567 bool result; 2568 2569 struct Intersection par; 2570 par.x1.CopyVector(&point1); 2571 par.x2.CopyVector(&point2); 2572 par.x3.CopyVector(&point3); 2573 par.x4.CopyVector(&point4); 2574 2575 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; 2576 gsl_multimin_fminimizer *s = NULL; 2577 gsl_vector *ss, *x; 2578 gsl_multimin_function minex_func; 2579 2580 size_t iter = 0; 2581 int status; 2582 double size; 2583 2584 /* Starting point */ 2585 x = gsl_vector_alloc(NDIM); 2586 gsl_vector_set(x, 0, point1.x[0]); 2587 gsl_vector_set(x, 1, point1.x[1]); 2588 gsl_vector_set(x, 2, point1.x[2]); 2589 2590 /* Set initial step sizes to 1 */ 2591 ss = gsl_vector_alloc(NDIM); 2592 gsl_vector_set_all(ss, 1.0); 2593 2594 /* Initialize method and iterate */ 2595 minex_func.n = NDIM; 2596 minex_func.f = &MinIntersectDistance; 2597 minex_func.params = (void *)∥ 2598 2599 s = gsl_multimin_fminimizer_alloc(T, NDIM); 2600 gsl_multimin_fminimizer_set(s, &minex_func, x, ss); 2601 2602 do { 2603 iter++; 2604 status = gsl_multimin_fminimizer_iterate(s); 2605 2606 if (status) { 2607 break; 2608 } 2609 2610 size = gsl_multimin_fminimizer_size(s); 2611 status = gsl_multimin_test_size(size, 1e-2); 2612 2613 if (status == GSL_SUCCESS) { 2614 cout << Verbose(2) << "converged to minimum" << endl; 2615 } 2616 } while (status == GSL_CONTINUE && iter < 100); 2617 2618 // check whether intersection is in between or not 2619 Vector intersection, SideA, SideB, HeightA, HeightB; 2620 double t1, t2; 2621 for (int i = 0; i < NDIM; i++) { 2622 intersection.x[i] = gsl_vector_get(s->x, i); 2623 } 2624 2625 SideA.CopyVector(&par.x2); 2626 SideA.SubtractVector(&par.x1); 2627 HeightA.CopyVector(&intersection); 2628 HeightA.SubtractVector(&par.x1); 2629 2630 t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA); 2631 2632 SideB.CopyVector(&par.x4); 2633 SideB.SubtractVector(&par.x3); 2634 HeightB.CopyVector(&intersection); 2635 HeightB.SubtractVector(&par.x3); 2636 2637 t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB); 2638 2639 cout << Verbose(2) << "Intersection " << intersection << " is at " 2640 << t1 << " for (" << point1 << "," << point2 << ") and at " 2641 << t2 << " for (" << point3 << "," << point4 << "): "; 2642 2643 if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) { 2644 cout << "true intersection." << endl; 2645 result = true; 2646 } else { 2647 cout << "intersection out of region of interest." << endl; 2648 result = false; 2649 } 2650 2651 // free minimizer stuff 2652 gsl_vector_free(x); 2653 gsl_vector_free(ss); 2654 gsl_multimin_fminimizer_free(s); 2655 2656 return result; 2657 } 2360 2658 2361 2659 /** Finds the second point of starting triangle. … … 2370 2668 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC) 2371 2669 { 2372 cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl; 2373 int i; 2374 Vector AngleCheck; 2375 atom* Walker; 2376 double norm = -1., angle; 2377 LinkedAtoms *List = NULL; 2378 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2379 2380 if (LC->SetIndexToAtom(a)) { // get cell for the starting atom 2381 for(int i=0;i<NDIM;i++) // store indices of this cell 2382 N[i] = LC->n[i]; 2383 } else { 2384 cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl; 2385 return; 2386 } 2387 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2388 cout << Verbose(2) << "LC Intervals:"; 2389 for (int i=0;i<NDIM;i++) { 2390 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2391 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2392 cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2393 } 2394 cout << endl; 2395 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2396 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2397 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2398 List = LC->GetCurrentCell(); 2399 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2400 if (List != NULL) { 2401 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2402 Candidate = (*Runner); 2403 // check if we only have one unique point yet ... 2404 if (a != Candidate) { 2405 cout << Verbose(3) << "Current candidate is " << *Candidate << ": "; 2406 AngleCheck.CopyVector(&(Candidate->x)); 2407 AngleCheck.SubtractVector(&(a->x)); 2408 norm = AngleCheck.Norm(); 2409 // second point shall have smallest angle with respect to Oben vector 2410 if (norm < RADIUS) { 2411 angle = AngleCheck.Angle(&Oben); 2412 if (angle < Storage[0]) { 2413 //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 2414 cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n"; 2415 Opt_Candidate = Candidate; 2416 Storage[0] = AngleCheck.Angle(&Oben); 2417 //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]); 2418 } else { 2419 cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl; 2420 } 2421 } else { 2422 cout << "Refused due to Radius " << norm << endl; 2423 } 2424 } 2425 } 2426 } 2427 } 2428 cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl; 2670 cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl; 2671 int i; 2672 Vector AngleCheck; 2673 atom* Walker; 2674 double norm = -1., angle; 2675 LinkedAtoms *List = NULL; 2676 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2677 2678 if (LC->SetIndexToAtom(a)) { // get cell for the starting atom 2679 for(int i=0;i<NDIM;i++) // store indices of this cell 2680 N[i] = LC->n[i]; 2681 } else { 2682 cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl; 2683 return; 2684 } 2685 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2686 cout << Verbose(2) << "LC Intervals from ["; 2687 for (int i=0;i<NDIM;i++) { 2688 cout << " " << N[i] << "<->" << LC->N[i]; 2689 } 2690 cout << "] :"; 2691 for (int i=0;i<NDIM;i++) { 2692 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2693 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2694 cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2695 } 2696 cout << endl; 2697 2698 2699 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2700 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2701 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2702 List = LC->GetCurrentCell(); 2703 cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2704 if (List != NULL) { 2705 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2706 Candidate = (*Runner); 2707 cout << Verbose(2) << "Current candidate is " << *Candidate << ": "; 2708 // check if we only have one unique point yet ... 2709 if (a != Candidate) { 2710 // Calculate center of the circle with radius RADIUS through points a and Candidate 2711 Vector OrthogonalizedOben, a_Candidate, Center; 2712 double distance, scaleFactor; 2713 2714 OrthogonalizedOben.CopyVector(&Oben); 2715 a_Candidate.CopyVector(&(a->x)); 2716 a_Candidate.SubtractVector(&(Candidate->x)); 2717 OrthogonalizedOben.ProjectOntoPlane(&a_Candidate); 2718 OrthogonalizedOben.Normalize(); 2719 distance = 0.5 * a_Candidate.Norm(); 2720 scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance))); 2721 OrthogonalizedOben.Scale(scaleFactor); 2722 2723 Center.CopyVector(&(Candidate->x)); 2724 Center.AddVector(&(a->x)); 2725 Center.Scale(0.5); 2726 Center.AddVector(&OrthogonalizedOben); 2727 2728 AngleCheck.CopyVector(&Center); 2729 AngleCheck.SubtractVector(&(a->x)); 2730 norm = a_Candidate.Norm(); 2731 // second point shall have smallest angle with respect to Oben vector 2732 if (norm < RADIUS) { 2733 angle = AngleCheck.Angle(&Oben); 2734 if (angle < Storage[0]) { 2735 //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 2736 cout << "Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n"; 2737 Opt_Candidate = Candidate; 2738 Storage[0] = angle; 2739 //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]); 2740 } else { 2741 cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl; 2742 } 2743 } else { 2744 cout << "Refused due to Radius " << norm << endl; 2745 } 2746 } else { 2747 cout << " Candidate is equal to first endpoint " << *a << "." << endl; 2748 } 2749 } 2750 } else { 2751 cout << "Linked cell list is empty." << endl; 2752 } 2753 } 2754 cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl; 2429 2755 }; 2430 2756 … … 2438 2764 void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC) 2439 2765 { 2440 cout << Verbose(1) << "Begin of Find_starting_triangle\n"; 2441 int i = 0; 2442 LinkedAtoms *List = NULL; 2443 atom* Walker; 2444 atom* FirstPoint; 2445 atom* SecondPoint; 2446 atom* MaxAtom[NDIM]; 2447 double max_coordinate[NDIM]; 2448 Vector Oben; 2449 Vector helper; 2450 Vector Chord; 2451 Vector SearchDirection; 2452 Vector OptCandidateCenter; 2453 2454 Oben.Zero(); 2455 2456 for (i = 0; i < 3; i++) { 2457 MaxAtom[i] = NULL; 2458 max_coordinate[i] = -1; 2459 } 2460 2461 // 1. searching topmost atom with respect to each axis 2462 for (int i=0;i<NDIM;i++) { // each axis 2463 LC->n[i] = LC->N[i]-1; // current axis is topmost cell 2464 for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++) 2465 for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) { 2466 List = LC->GetCurrentCell(); 2467 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2468 if (List != NULL) { 2469 for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) { 2470 cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl; 2471 if ((*Runner)->x.x[i] > max_coordinate[i]) { 2472 max_coordinate[i] = (*Runner)->x.x[i]; 2473 MaxAtom[i] = (*Runner); 2474 } 2475 } 2476 } else { 2477 cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl; 2478 } 2479 } 2480 } 2481 2482 cout << Verbose(2) << "Found maximum coordinates: "; 2483 for (int i=0;i<NDIM;i++) 2484 cout << i << ": " << *MaxAtom[i] << "\t"; 2485 cout << endl; 2486 const int k = 1; // arbitrary choice 2487 Oben.x[k] = 1.; 2488 FirstPoint = MaxAtom[k]; 2489 cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl; 2490 2491 // add first point 2492 AddTrianglePoint(FirstPoint, 0); 2493 2494 double ShortestAngle; 2495 atom* Opt_Candidate = NULL; 2496 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant. 2497 2498 Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 2499 SecondPoint = Opt_Candidate; 2500 cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n"; 2501 2502 // add second point and first baseline 2503 AddTrianglePoint(SecondPoint, 1); 2504 AddTriangleLine(TPS[0], TPS[1], 0); 2505 2506 helper.CopyVector(&(FirstPoint->x)); 2507 helper.SubtractVector(&(SecondPoint->x)); 2508 helper.Normalize(); 2509 Oben.ProjectOntoPlane(&helper); 2510 Oben.Normalize(); 2511 helper.VectorProduct(&Oben); 2512 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2513 2514 Chord.CopyVector(&(FirstPoint->x)); // bring into calling function 2515 Chord.SubtractVector(&(SecondPoint->x)); 2516 double radius = Chord.ScalarProduct(&Chord); 2517 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 2518 helper.CopyVector(&Oben); 2519 helper.Scale(CircleRadius); 2520 // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized) 2521 2522 cout << Verbose(2) << "Looking for third point candidates \n"; 2523 // look in one direction of baseline for initial candidate 2524 Opt_Candidate = NULL; 2525 SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ... 2526 2527 cout << Verbose(1) << "Looking for third point candidates ...\n"; 2528 Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC); 2529 cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl; 2530 2531 // add third point 2532 AddTrianglePoint(Opt_Candidate, 2); 2533 2534 // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate 2535 2536 // Finally, we only have to add the found further lines 2537 AddTriangleLine(TPS[1], TPS[2], 1); 2538 AddTriangleLine(TPS[0], TPS[2], 2); 2539 // ... and triangles to the Maps of the Tesselation class 2540 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2541 AddTriangleToLines(); 2542 // ... and calculate its normal vector (with correct orientation) 2543 OptCandidateCenter.Scale(-1.); 2544 cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl; 2545 BTS->GetNormalVector(OptCandidateCenter); 2546 cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n"; 2547 cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl; 2548 cout << Verbose(1) << "End of Find_starting_triangle\n"; 2766 cout << Verbose(1) << "Begin of Find_starting_triangle\n"; 2767 int i = 0; 2768 LinkedAtoms *List = NULL; 2769 atom* Walker; 2770 atom* FirstPoint; 2771 atom* SecondPoint; 2772 atom* MaxAtom[NDIM]; 2773 double max_coordinate[NDIM]; 2774 Vector Oben; 2775 Vector helper; 2776 Vector Chord; 2777 Vector SearchDirection; 2778 2779 Oben.Zero(); 2780 2781 for (i = 0; i < 3; i++) { 2782 MaxAtom[i] = NULL; 2783 max_coordinate[i] = -1; 2784 } 2785 2786 // 1. searching topmost atom with respect to each axis 2787 for (int i=0;i<NDIM;i++) { // each axis 2788 LC->n[i] = LC->N[i]-1; // current axis is topmost cell 2789 for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++) 2790 for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) { 2791 List = LC->GetCurrentCell(); 2792 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2793 if (List != NULL) { 2794 for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) { 2795 cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl; 2796 if ((*Runner)->x.x[i] > max_coordinate[i]) { 2797 max_coordinate[i] = (*Runner)->x.x[i]; 2798 MaxAtom[i] = (*Runner); 2799 } 2800 } 2801 } else { 2802 cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl; 2803 } 2804 } 2805 } 2806 2807 cout << Verbose(2) << "Found maximum coordinates: "; 2808 for (int i=0;i<NDIM;i++) 2809 cout << i << ": " << *MaxAtom[i] << "\t"; 2810 cout << endl; 2811 const int k = 1; // arbitrary choice 2812 Oben.x[k] = 1.; 2813 FirstPoint = MaxAtom[k]; 2814 cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl; 2815 2816 double ShortestAngle; 2817 atom* Opt_Candidate = NULL; 2818 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant. 2819 2820 Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 2821 SecondPoint = Opt_Candidate; 2822 cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n"; 2823 2824 helper.CopyVector(&(FirstPoint->x)); 2825 helper.SubtractVector(&(SecondPoint->x)); 2826 helper.Normalize(); 2827 Oben.ProjectOntoPlane(&helper); 2828 Oben.Normalize(); 2829 helper.VectorProduct(&Oben); 2830 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2831 2832 Chord.CopyVector(&(FirstPoint->x)); // bring into calling function 2833 Chord.SubtractVector(&(SecondPoint->x)); 2834 double radius = Chord.ScalarProduct(&Chord); 2835 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 2836 helper.CopyVector(&Oben); 2837 helper.Scale(CircleRadius); 2838 // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized) 2839 2840 cout << Verbose(2) << "Looking for third point candidates \n"; 2841 // look in one direction of baseline for initial candidate 2842 CandidateList *Opt_Candidates = new CandidateList(); 2843 SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ... 2844 2845 // adding point 1 and point 2 and the line between them 2846 AddTrianglePoint(FirstPoint, 0); 2847 AddTrianglePoint(SecondPoint, 1); 2848 AddTriangleLine(TPS[0], TPS[1], 0); 2849 2850 cout << Verbose(1) << "Looking for third point candidates ...\n"; 2851 cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n"; 2852 Find_third_point_for_Tesselation( 2853 Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC 2854 ); 2855 cout << Verbose(1) << "Third Points are "; 2856 CandidateList::iterator it; 2857 for (it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2858 cout << " " << *(*it)->point; 2859 } 2860 cout << endl; 2861 2862 for (it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2863 // add third triangle point 2864 AddTrianglePoint((*it)->point, 2); 2865 // add the second and third line 2866 AddTriangleLine(TPS[1], TPS[2], 1); 2867 AddTriangleLine(TPS[0], TPS[2], 2); 2868 // ... and triangles to the Maps of the Tesselation class 2869 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2870 AddTriangleToLines(); 2871 // ... and calculate its normal vector (with correct orientation) 2872 (*it)->OptCenter.Scale(-1.); 2873 cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl; 2874 BTS->GetNormalVector((*it)->OptCenter); // vector to compare with should point inwards 2875 cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " 2876 << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n"; 2877 2878 // if we do not reach the end with the next step of iteration, we need to setup a new first line 2879 if (it != Opt_Candidates->end()--) { 2880 FirstPoint = (*it)->BaseLine->endpoints[0]->node; 2881 SecondPoint = (*it)->point; 2882 // adding point 1 and point 2 and the line between them 2883 AddTrianglePoint(FirstPoint, 0); 2884 AddTrianglePoint(SecondPoint, 1); 2885 AddTriangleLine(TPS[0], TPS[1], 0); 2886 } 2887 } 2888 cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl; 2889 cout << Verbose(1) << "End of Find_starting_triangle\n"; 2549 2890 }; 2550 2891 … … 2560 2901 */ 2561 2902 bool Tesselation::Find_next_suitable_triangle(ofstream *out, 2562 molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, 2563 const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC) 2564 { 2565 cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n"; 2566 ofstream *tempstream = NULL; 2567 char NumberName[255]; 2568 double tmp; 2569 2570 atom* Opt_Candidate = NULL; 2571 Vector OptCandidateCenter; 2572 2573 Vector CircleCenter; 2574 Vector CirclePlaneNormal; 2575 Vector OldSphereCenter; 2576 Vector SearchDirection; 2577 Vector helper; 2578 atom *ThirdNode = NULL; 2579 double ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2580 double radius, CircleRadius; 2581 2582 cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl; 2583 for (int i=0;i<3;i++) 2584 if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node)) 2585 ThirdNode = T.endpoints[i]->node; 2586 2587 // construct center of circle 2588 CircleCenter.CopyVector(&Line.endpoints[0]->node->x); 2589 CircleCenter.AddVector(&Line.endpoints[1]->node->x); 2590 CircleCenter.Scale(0.5); 2591 2592 // construct normal vector of circle 2593 CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x); 2594 CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x); 2595 2596 // calculate squared radius of circle 2597 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2598 if (radius/4. < RADIUS*RADIUS) { 2599 CircleRadius = RADIUS*RADIUS - radius/4.; 2600 CirclePlaneNormal.Normalize(); 2601 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2602 2603 // construct old center 2604 GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x)); 2605 helper.CopyVector(&T.NormalVector); // normal vector ensures that this is correct center of the two possible ones 2606 radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter); 2607 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2608 OldSphereCenter.AddVector(&helper); 2609 OldSphereCenter.SubtractVector(&CircleCenter); 2610 cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2611 2612 // construct SearchDirection 2613 SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal); 2614 helper.CopyVector(&Line.endpoints[0]->node->x); 2615 helper.SubtractVector(&ThirdNode->x); 2616 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards! 2617 SearchDirection.Scale(-1.); 2618 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2619 SearchDirection.Normalize(); 2620 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2621 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2622 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2623 } 2624 2625 // add third point 2626 cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl; 2627 Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC); 2628 2629 } else { 2630 cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl; 2903 molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, 2904 const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC) 2905 { 2906 cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n"; 2907 ofstream *tempstream = NULL; 2908 char NumberName[255]; 2909 double tmp; 2910 bool result = true; 2911 CandidateList *Opt_Candidates = new CandidateList(); 2912 2913 Vector CircleCenter; 2914 Vector CirclePlaneNormal; 2915 Vector OldSphereCenter; 2916 Vector SearchDirection; 2917 Vector helper; 2918 atom *ThirdNode = NULL; 2919 LineMap::iterator testline; 2920 double ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2921 double radius, CircleRadius; 2922 2923 cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl; 2924 for (int i=0;i<3;i++) 2925 if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node)) 2926 ThirdNode = T.endpoints[i]->node; 2927 2928 // construct center of circle 2929 CircleCenter.CopyVector(&Line.endpoints[0]->node->x); 2930 CircleCenter.AddVector(&Line.endpoints[1]->node->x); 2931 CircleCenter.Scale(0.5); 2932 2933 // construct normal vector of circle 2934 CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x); 2935 CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x); 2936 2937 // calculate squared radius of circle 2938 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2939 if (radius/4. < RADIUS*RADIUS) { 2940 CircleRadius = RADIUS*RADIUS - radius/4.; 2941 CirclePlaneNormal.Normalize(); 2942 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2943 2944 // construct old center 2945 GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x)); 2946 helper.CopyVector(&T.NormalVector); // normal vector ensures that this is correct center of the two possible ones 2947 radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter); 2948 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2949 OldSphereCenter.AddVector(&helper); 2950 OldSphereCenter.SubtractVector(&CircleCenter); 2951 cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2952 2953 // construct SearchDirection 2954 SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal); 2955 helper.CopyVector(&Line.endpoints[0]->node->x); 2956 helper.SubtractVector(&ThirdNode->x); 2957 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 2958 SearchDirection.Scale(-1.); 2959 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2960 SearchDirection.Normalize(); 2961 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2962 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2963 // rotated the wrong way! 2964 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2965 } 2966 2967 // add third point 2968 cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl; 2969 Find_third_point_for_Tesselation( 2970 T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates, 2971 &ShortestAngle, RADIUS, LC 2972 ); 2973 2974 } else { 2975 cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl; 2976 } 2977 2978 if (Opt_Candidates->begin() == Opt_Candidates->end()) { 2979 cerr << "WARNING: Could not find a suitable candidate." << endl; 2980 return false; 2981 } 2982 cout << Verbose(1) << "Third Points are "; 2983 CandidateList::iterator it; 2984 for (it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2985 cout << " " << *(*it)->point; 2986 } 2987 cout << endl; 2988 2989 BoundaryLineSet *BaseRay = &Line; 2990 for (it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2991 cout << Verbose(1) << " Third point candidate is " << *(*it)->point 2992 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2993 cout << Verbose(1) << " Baseline is " << BaseRay << endl; 2994 2995 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2996 atom *AtomCandidates[3]; 2997 AtomCandidates[0] = (*it)->point; 2998 AtomCandidates[1] = BaseRay->endpoints[0]->node; 2999 AtomCandidates[2] = BaseRay->endpoints[1]->node; 3000 int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates); 3001 3002 BTS = NULL; 3003 // If there is no triangle, add it regularly. 3004 if (existentTrianglesCount == 0) { 3005 AddTrianglePoint((*it)->point, 0); 3006 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); 3007 AddTrianglePoint(BaseRay->endpoints[1]->node, 2); 3008 3009 AddTriangleLine(TPS[0], TPS[1], 0); 3010 AddTriangleLine(TPS[0], TPS[2], 1); 3011 AddTriangleLine(TPS[1], TPS[2], 2); 3012 3013 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 3014 AddTriangleToLines(); 3015 (*it)->OptCenter.Scale(-1.); 3016 BTS->GetNormalVector((*it)->OptCenter); 3017 (*it)->OptCenter.Scale(-1.); 3018 3019 cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector 3020 << " for this triangle ... " << endl; 3021 cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << BaseRay << "." << endl; 3022 } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time. 3023 AddTrianglePoint((*it)->point, 0); 3024 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); 3025 AddTrianglePoint(BaseRay->endpoints[1]->node, 2); 3026 3027 AddTriangleLine(TPS[0], TPS[1], 0); 3028 AddTriangleLine(TPS[0], TPS[2], 1); 3029 AddTriangleLine(TPS[1], TPS[2], 2); 3030 3031 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 3032 //TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 3033 AddTriangleToLines(); 3034 3035 (*it)->OtherOptCenter.Scale(-1.); 3036 BTS->GetNormalVector((*it)->OtherOptCenter); 3037 (*it)->OtherOptCenter.Scale(-1.); 3038 3039 cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector 3040 << " for this triangle ... " << endl; 3041 cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl; 3042 } else { 3043 cout << Verbose(1) << "This triangle consisting of "; 3044 cout << *(*it)->point << ", "; 3045 cout << *BaseRay->endpoints[0]->node << " and "; 3046 cout << *BaseRay->endpoints[1]->node << " "; 3047 cout << "is invalid!" << endl; 3048 result = false; 3049 } 3050 3051 if ((existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration 3052 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name); 3053 if (DoTecplotOutput) { 3054 string NameofTempFile(tempbasename); 3055 NameofTempFile.append(NumberName); 3056 for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos)) 3057 NameofTempFile.erase(npos, 1); 3058 NameofTempFile.append(TecplotSuffix); 3059 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 3060 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 3061 write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten); 3062 tempstream->close(); 3063 tempstream->flush(); 3064 delete(tempstream); 3065 } 3066 3067 if (DoRaster3DOutput) { 3068 string NameofTempFile(tempbasename); 3069 NameofTempFile.append(NumberName); 3070 for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos)) 3071 NameofTempFile.erase(npos, 1); 3072 NameofTempFile.append(Raster3DSuffix); 3073 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 3074 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 3075 write_raster3d_file(out, tempstream, this, mol); 3076 // include the current position of the virtual sphere in the temporary raster3d file 3077 // make the circumsphere's center absolute again 3078 helper.CopyVector(&BaseRay->endpoints[0]->node->x); 3079 helper.AddVector(&BaseRay->endpoints[1]->node->x); 3080 helper.Scale(0.5); 3081 (*it)->OptCenter.AddVector(&helper); 3082 Vector *center = mol->DetermineCenterOfAll(out); 3083 (*it)->OptCenter.AddVector(center); 3084 delete(center); 3085 // and add to file plus translucency object 3086 *tempstream << "# current virtual sphere\n"; 3087 *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 3088 *tempstream << "2\n " << (*it)->OptCenter.x[0] << " " 3089 << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2] 3090 << "\t" << RADIUS << "\t1 0 0\n"; 3091 *tempstream << "9\n terminating special property\n"; 3092 tempstream->close(); 3093 tempstream->flush(); 3094 delete(tempstream); 3095 } 3096 if (DoTecplotOutput || DoRaster3DOutput) 3097 TriangleFilesWritten++; 3098 } 3099 3100 // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 3101 BaseRay = BLS[0]; 3102 // LineMap::iterator LineIterator = Line.endpoints[0]->lines.find((*it)->point->nr); 3103 // for (; LineIterator != Line.endpoints[0]->lines.end(); LineIterator++) { 3104 // if ((*LineIterator->second).TrianglesCount != 2) 3105 // break; 3106 // } 3107 // if (LineIterator == Line.endpoints[0]->lines.end()) 3108 // cout << Verbose(1) << "ERROR: I could not find a suitable line with less than two triangles connected!" << endl; 3109 } 3110 3111 cout << Verbose(1) << "End of Find_next_suitable_triangle\n"; 3112 return result; 3113 }; 3114 3115 /** 3116 * Sort function for the candidate list. 3117 */ 3118 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) { 3119 Vector BaseLineVector, OrthogonalVector, helper; 3120 if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 3121 cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 3122 //return false; 3123 exit(1); 2631 3124 } 2632 2633 if (Opt_Candidate == NULL) { 2634 cerr << "WARNING: Could not find a suitable candidate." << endl; 2635 return false; 3125 // create baseline vector 3126 BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x)); 3127 BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 3128 BaseLineVector.Normalize(); 3129 3130 // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 3131 helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 3132 helper.SubtractVector(&(candidate1->point->x)); 3133 OrthogonalVector.CopyVector(&helper); 3134 helper.VectorProduct(&BaseLineVector); 3135 OrthogonalVector.SubtractVector(&helper); 3136 OrthogonalVector.Normalize(); 3137 3138 // calculate both angles and correct with in-plane vector 3139 helper.CopyVector(&(candidate1->point->x)); 3140 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 3141 double phi = BaseLineVector.Angle(&helper); 3142 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 3143 phi = 2.*M_PI - phi; 3144 } 3145 helper.CopyVector(&(candidate2->point->x)); 3146 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 3147 double psi = BaseLineVector.Angle(&helper); 3148 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 3149 psi = 2.*M_PI - psi; 2636 3150 } 2637 cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl; 2638 2639 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2640 atom *AtomCandidates[3]; 2641 AtomCandidates[0] = Opt_Candidate; 2642 AtomCandidates[1] = Line.endpoints[0]->node; 2643 AtomCandidates[2] = Line.endpoints[1]->node; 2644 bool flag = CheckPresenceOfTriangle(out, AtomCandidates); 2645 2646 if (flag) { // if so, add 2647 AddTrianglePoint(Opt_Candidate, 0); 2648 AddTrianglePoint(Line.endpoints[0]->node, 1); 2649 AddTrianglePoint(Line.endpoints[1]->node, 2); 2650 2651 AddTriangleLine(TPS[0], TPS[1], 0); 2652 AddTriangleLine(TPS[0], TPS[2], 1); 2653 AddTriangleLine(TPS[1], TPS[2], 2); 2654 2655 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2656 AddTriangleToLines(); 2657 2658 OptCandidateCenter.Scale(-1.); 2659 BTS->GetNormalVector(OptCandidateCenter); 2660 OptCandidateCenter.Scale(-1.); 2661 2662 cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl; 2663 cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl; 2664 } else { // else, yell and do nothing 2665 cout << Verbose(1) << "This triangle consisting of "; 2666 cout << *Opt_Candidate << ", "; 2667 cout << *Line.endpoints[0]->node << " and "; 2668 cout << *Line.endpoints[1]->node << " "; 2669 cout << "is invalid!" << endl; 2670 return false; 2671 } 2672 2673 if (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration 2674 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name); 2675 if (DoTecplotOutput) { 2676 string NameofTempFile(tempbasename); 2677 NameofTempFile.append(NumberName); 2678 for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos)) 2679 NameofTempFile.erase(npos, 1); 2680 NameofTempFile.append(TecplotSuffix); 2681 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 2682 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 2683 write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten); 2684 tempstream->close(); 2685 tempstream->flush(); 2686 delete(tempstream); 2687 } 2688 2689 if (DoRaster3DOutput) { 2690 string NameofTempFile(tempbasename); 2691 NameofTempFile.append(NumberName); 2692 for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos)) 2693 NameofTempFile.erase(npos, 1); 2694 NameofTempFile.append(Raster3DSuffix); 2695 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 2696 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 2697 write_raster3d_file(out, tempstream, this, mol); 2698 // include the current position of the virtual sphere in the temporary raster3d file 2699 // make the circumsphere's center absolute again 2700 helper.CopyVector(&Line.endpoints[0]->node->x); 2701 helper.AddVector(&Line.endpoints[1]->node->x); 2702 helper.Scale(0.5); 2703 OptCandidateCenter.AddVector(&helper); 2704 Vector *center = mol->DetermineCenterOfAll(out); 2705 OptCandidateCenter.AddVector(center); 2706 delete(center); 2707 // and add to file plus translucency object 2708 *tempstream << "# current virtual sphere\n"; 2709 *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 2710 *tempstream << "2\n " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n"; 2711 *tempstream << "9\n terminating special property\n"; 2712 tempstream->close(); 2713 tempstream->flush(); 2714 delete(tempstream); 2715 } 2716 if (DoTecplotOutput || DoRaster3DOutput) 2717 TriangleFilesWritten++; 2718 } 2719 2720 cout << Verbose(1) << "End of Find_next_suitable_triangle\n"; 2721 return true; 2722 }; 3151 3152 cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl; 3153 cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl; 3154 3155 // return comparison 3156 return phi < psi; 3157 } 2723 3158 2724 3159 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule. … … 2726 3161 * \param *mol molecule structure with Atom's and Bond's 2727 3162 * \param *Tess Tesselation filled with points, lines and triangles on boundary on return 2728 * \param *LCList linked cell list of all atoms2729 3163 * \param *filename filename prefix for output of vertex data 2730 3164 * \para RADIUS radius of the virtual sphere … … 2732 3166 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS) 2733 3167 { 2734 int N = 0; 2735 bool freeTess = false; 2736 *out << Verbose(1) << "Entering search for non convex hull. " << endl; 2737 if (Tess == NULL) { 2738 *out << Verbose(1) << "Allocating Tesselation struct ..." << endl; 2739 Tess = new Tesselation; 2740 freeTess = true; 2741 } 2742 bool freeLC = false; 2743 LineMap::iterator baseline; 2744 *out << Verbose(0) << "Begin of Find_non_convex_border\n"; 2745 bool flag = false; // marks whether we went once through all baselines without finding any without two triangles 2746 bool failflag = false; 2747 2748 if (LCList == NULL) { 2749 LCList = new LinkedCell(mol, 2.*RADIUS); 2750 freeLC = true; 2751 } 2752 2753 Tess->Find_starting_triangle(out, mol, RADIUS, LCList); 2754 2755 baseline = Tess->LinesOnBoundary.begin(); 2756 while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) { 2757 if (baseline->second->TrianglesCount == 1) { 2758 failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one. 2759 flag = flag || failflag; 2760 if (!failflag) 2761 cerr << "WARNING: Find_next_suitable_triangle failed." << endl; 2762 } else { 2763 cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl; 2764 } 2765 N++; 2766 baseline++; 2767 if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) { 2768 baseline = Tess->LinesOnBoundary.begin(); // restart if we reach end due to newly inserted lines 2769 flag = false; 2770 } 2771 } 2772 if (1) { //failflag) { 2773 *out << Verbose(1) << "Writing final tecplot file\n"; 2774 if (DoTecplotOutput) { 2775 string OutputName(filename); 2776 OutputName.append(TecplotSuffix); 2777 ofstream *tecplot = new ofstream(OutputName.c_str()); 2778 write_tecplot_file(out, tecplot, Tess, mol, -1); 2779 tecplot->close(); 2780 delete(tecplot); 2781 } 2782 if (DoRaster3DOutput) { 2783 string OutputName(filename); 2784 OutputName.append(Raster3DSuffix); 2785 ofstream *raster = new ofstream(OutputName.c_str()); 2786 write_raster3d_file(out, raster, Tess, mol); 2787 raster->close(); 2788 delete(raster); 2789 } 2790 } else { 2791 cerr << "ERROR: Could definately not find all necessary triangles!" << endl; 2792 } 2793 if (freeTess) 2794 delete(Tess); 2795 if (freeLC) 2796 delete(LCList); 2797 *out << Verbose(0) << "End of Find_non_convex_border\n"; 3168 int N = 0; 3169 bool freeTess = false; 3170 bool freeLC = false; 3171 *out << Verbose(1) << "Entering search for non convex hull. " << endl; 3172 if (Tess == NULL) { 3173 *out << Verbose(1) << "Allocating Tesselation struct ..." << endl; 3174 Tess = new Tesselation; 3175 freeTess = true; 3176 } 3177 LineMap::iterator baseline; 3178 LineMap::iterator testline; 3179 *out << Verbose(0) << "Begin of Find_non_convex_border\n"; 3180 bool flag = false; // marks whether we went once through all baselines without finding any without two triangles 3181 bool failflag = false; 3182 3183 if (LCList == NULL) { 3184 LCList = new LinkedCell(mol, 2.*RADIUS); 3185 freeLC = true; 3186 } 3187 3188 Tess->Find_starting_triangle(out, mol, RADIUS, LCList); 3189 3190 baseline = Tess->LinesOnBoundary.begin(); 3191 while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) { 3192 if (baseline->second->TrianglesCount == 1) { 3193 failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one. 3194 flag = flag || failflag; 3195 if (!failflag) 3196 cerr << "WARNING: Find_next_suitable_triangle failed." << endl; 3197 3198 // we inserted new lines, hence show list with connected triangles 3199 cout << Verbose(1) << "List of Baselines with connected triangles so far:" << endl; 3200 for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) { 3201 cout << Verbose(1) << *testline->second << "\t" << testline->second->TrianglesCount << endl; 3202 } 3203 } else { 3204 cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl; 3205 if (baseline->second->TrianglesCount != 2) 3206 cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl; 3207 } 3208 3209 N++; 3210 baseline++; 3211 if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) { 3212 baseline = Tess->LinesOnBoundary.begin(); // restart if we reach end due to newly inserted lines 3213 flag = false; 3214 } 3215 } 3216 if (1) { //failflag) { 3217 *out << Verbose(1) << "Writing final tecplot file\n"; 3218 if (DoTecplotOutput) { 3219 string OutputName(filename); 3220 OutputName.append(TecplotSuffix); 3221 ofstream *tecplot = new ofstream(OutputName.c_str()); 3222 write_tecplot_file(out, tecplot, Tess, mol, -1); 3223 tecplot->close(); 3224 delete(tecplot); 3225 } 3226 if (DoRaster3DOutput) { 3227 string OutputName(filename); 3228 OutputName.append(Raster3DSuffix); 3229 ofstream *raster = new ofstream(OutputName.c_str()); 3230 write_raster3d_file(out, raster, Tess, mol); 3231 raster->close(); 3232 delete(raster); 3233 } 3234 } else { 3235 cerr << "ERROR: Could definitively not find all necessary triangles!" << endl; 3236 } 3237 if (freeTess) 3238 delete(Tess); 3239 if (freeLC) 3240 delete(LCList); 3241 *out << Verbose(0) << "End of Find_non_convex_border\n"; 2798 3242 }; 2799 3243 -
src/boundary.hpp
rca2587 r3d919e 5 5 class BoundaryLineSet; 6 6 class BoundaryTriangleSet; 7 class CandidateForTesselation; 7 8 8 9 // include config.h … … 17 18 18 19 #include <gsl/gsl_poly.h> 19 #include <gsl/gsl_matrix.h>20 #include <gsl/gsl_linalg.h>21 #include <gsl/gsl_multimin.h>22 #include <gsl/gsl_permutation.h>23 20 24 21 #include "linkedcell.hpp" … … 27 24 template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2) 28 25 { 29 30 31 32 33 34 35 26 if (endpoint1->Nr < endpoint2->Nr) { 27 endpoints[0] = endpoint1; 28 endpoints[1] = endpoint2; 29 } else { 30 endpoints[0] = endpoint2; 31 endpoints[1] = endpoint1; 32 } 36 33 }; 37 34 38 35 class BoundaryPointSet { 39 40 41 42 36 public: 37 BoundaryPointSet(); 38 BoundaryPointSet(atom *Walker); 39 ~BoundaryPointSet(); 43 40 44 41 void AddLine(class BoundaryLineSet *line); 45 42 46 47 48 49 43 LineMap lines; 44 int LinesCount; 45 atom *node; 46 int Nr; 50 47 }; 51 48 52 49 class BoundaryLineSet { 53 54 55 56 50 public: 51 BoundaryLineSet(); 52 BoundaryLineSet(class BoundaryPointSet *Point[2], int number); 53 ~BoundaryLineSet(); 57 54 58 55 void AddTriangle(class BoundaryTriangleSet *triangle); 59 56 60 61 62 63 57 class BoundaryPointSet *endpoints[2]; 58 TriangleMap triangles; 59 int TrianglesCount; 60 int Nr; 64 61 }; 65 62 66 63 class BoundaryTriangleSet { 67 68 69 70 64 public: 65 BoundaryTriangleSet(); 66 BoundaryTriangleSet(class BoundaryLineSet *line[3], int number); 67 ~BoundaryTriangleSet(); 71 68 72 69 void GetNormalVector(Vector &NormalVector); 73 70 74 75 76 77 71 class BoundaryPointSet *endpoints[3]; 72 class BoundaryLineSet *lines[3]; 73 Vector NormalVector; 74 int Nr; 78 75 }; 79 76 77 78 class CandidateForTesselation { 79 public : 80 CandidateForTesselation(atom* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter); 81 ~CandidateForTesselation(); 82 atom *point; 83 BoundaryLineSet *BaseLine; 84 Vector OptCenter; 85 Vector OtherOptCenter; 86 }; 87 88 80 89 class Tesselation { 81 90 public: 82 91 83 84 92 Tesselation(); 93 ~Tesselation(); 85 94 86 void TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol); 87 void GuessStartingTriangle(ofstream *out); 88 void AddPoint(atom * Walker); 89 void AddTrianglePoint(atom* Candidate, int n); 90 void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n); 91 void AddTriangleToLines(); 92 void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC); 93 bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC); 94 bool CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]); 95 void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol); 95 void TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol); 96 void GuessStartingTriangle(ofstream *out); 97 void AddPoint(atom * Walker); 98 void AddTrianglePoint(atom* Candidate, int n); 99 void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n); 100 void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n); 101 void AddTriangleToLines(); 102 void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC); 103 bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC); 104 int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]); 105 void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol); 96 106 97 98 99 100 101 102 103 104 105 106 107 107 PointMap PointsOnBoundary; 108 LineMap LinesOnBoundary; 109 TriangleMap TrianglesOnBoundary; 110 class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions 111 class BoundaryPointSet *BPS[2]; 112 class BoundaryLineSet *BLS[3]; 113 class BoundaryTriangleSet *BTS; 114 int PointsOnBoundaryCount; 115 int LinesOnBoundaryCount; 116 int TrianglesOnBoundaryCount; 117 int TriangleFilesWritten; 108 118 }; 109 119 … … 117 127 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem); 118 128 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity); 119 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC List, const char *tempbasename, const double RADIUS);129 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS); 120 130 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC); 121 131 bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess); 132 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4); 133 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2); 122 134 123 135 #endif /*BOUNDARY_HPP_*/ -
src/molecules.hpp
rca2587 r3d919e 57 57 #define PointPair pair < int, class BoundaryPointSet * > 58 58 #define PointTestPair pair < PointMap::iterator, bool > 59 60 #define LineMap map < int, class BoundaryLineSet * > 59 #define CandidateList list <class CandidateForTesselation *> 60 61 #define LineMap multimap < int, class BoundaryLineSet * > 61 62 #define LinePair pair < int, class BoundaryLineSet * > 62 63 #define LineTestPair pair < LineMap::iterator, bool > … … 71 72 #define MoleculeList list <molecule *> 72 73 #define MoleculeListTest pair <MoleculeList::iterator, bool> 74 75 #define LinkedAtoms list <atom *> 73 76 74 77 /******************************** Some small functions and/or structures **********************************/
Note:
See TracChangeset
for help on using the changeset viewer.