- Timestamp:
- Jul 7, 2009, 6:49:32 AM (16 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, Candidate_v1.7.0, 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
-
boundary.cpp (modified) (46 diffs)
-
boundary.hpp (modified) (4 diffs)
-
molecules.hpp (modified) (2 diffs)
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 LinesCount = 0;18 Nr = -1;23 LinesCount = 0; 24 Nr = -1; 19 25 } 20 26 ; … … 22 28 BoundaryPointSet::BoundaryPointSet(atom *Walker) 23 29 { 24 node = Walker;25 LinesCount = 0;26 Nr = Walker->nr;30 node = Walker; 31 LinesCount = 0; 32 Nr = Walker->nr; 27 33 } 28 34 ; … … 30 36 BoundaryPointSet::~BoundaryPointSet() 31 37 { 32 cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;33 if (!lines.empty())34 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;35 node = NULL;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 cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."42 << endl;43 if (line->endpoints[0] == this)44 {45 lines.insert(LinePair(line->endpoints[1]->Nr, line));46 }47 else48 {49 lines.insert(LinePair(line->endpoints[0]->Nr, line));50 }51 LinesCount++;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 ost << "[" << a.Nr << "|" << a.node->Name << "]";59 return ost;64 ost << "[" << a.Nr << "|" << a.node->Name << "]"; 65 return ost; 60 66 } 61 67 ; … … 65 71 BoundaryLineSet::BoundaryLineSet() 66 72 { 67 for (int i = 0; i < 2; i++)68 endpoints[i] = NULL;69 TrianglesCount = 0;70 Nr = -1;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 // set number77 Nr = number;78 // set endpoints in ascending order79 SetEndpointsOrdered(endpoints, Point[0], Point[1]);80 // add this line to the hash maps of both endpoints81 Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.82 Point[1]->AddLine(this); //83 // clear triangles list84 TrianglesCount = 0;85 cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;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 int Numbers[2];92 Numbers[0] = endpoints[1]->Nr;93 Numbers[1] = endpoints[0]->Nr;94 for (int i = 0; i < 2; i++) {95 cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;96 endpoints[i]->lines.erase(Numbers[i]);97 if (endpoints[i]->lines.empty()) {98 cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;99 if (endpoints[i] != NULL) {100 delete(endpoints[i]);101 endpoints[i] = NULL;102 } else103 cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;104 } else105 cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;106 }107 if (!triangles.empty())108 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;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 cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."116 << endl;117 triangles.insert(TrianglePair(triangle->Nr, triangle));118 TrianglesCount++;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 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","126 << a.endpoints[1]->node->Name << "]";127 return ost;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 for (int i = 0; i < 3; i++)137 {138 endpoints[i] = NULL;139 lines[i] = NULL;140 }141 Nr = -1;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 for (int i = 0; i < 3; i++) {190 cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;191 lines[i]->triangles.erase(Nr);192 if (lines[i]->triangles.empty()) {193 cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;194 if (lines[i] != NULL) {195 delete (lines[i]);196 lines[i] = NULL;197 } else198 cerr << "ERROR: This line " << i << " has already been free'd." << endl;199 } else200 cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;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 // get normal vector209 NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,210 &endpoints[2]->node->x);211 212 // make it always point inward (any offset vector onto plane projected onto normal vector suffices)213 if (NormalVector.Projection(&OtherVector) > 0)214 NormalVector.Scale(-1.);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 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","222 << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";223 return ost;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 class BoundaryLineSet * lines[2] =238 { line1, line2 };239 class BoundaryPointSet *node = NULL;240 map<int, class BoundaryPointSet *> OrderMap;241 pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;242 for (int i = 0; i < 2; i++)243 // for both lines244 for (int j = 0; j < 2; j++)245 { // for both endpoints246 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (247 lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));248 if (!OrderTest.second)249 { // if insertion fails, we have common endpoint250 node = OrderTest.first->second;251 cout << Verbose(5) << "Common endpoint of lines " << *line1252 << " and " << *line2 << " is: " << *node << "." << endl;253 j = 2;254 i = 2;255 break;256 }257 }258 return node;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 atom *Walker = NULL;273 PointMap PointsOnBoundary;274 LineMap LinesOnBoundary;275 TriangleMap TrianglesOnBoundary;276 277 *out << Verbose(1) << "Finding all boundary points." << endl;278 Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)279 BoundariesTestPair BoundaryTestPair;280 Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;281 double radius, angle;282 // 3a. Go through every axis283 for (int axis = 0; axis < NDIM; axis++)284 {285 AxisVector.Zero();286 AngleReferenceVector.Zero();287 AngleReferenceNormalVector.Zero();288 AxisVector.x[axis] = 1.;289 AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;290 AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;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 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours300 Walker = mol->start;301 while (Walker->next != mol->end)302 {303 Walker = Walker->next;304 Vector ProjectedVector;305 ProjectedVector.CopyVector(&Walker->x);306 ProjectedVector.ProjectOntoPlane(&AxisVector);307 // correct for negative side308 //if (Projection(y) < 0)309 //angle = 2.*M_PI - angle;310 radius = ProjectedVector.Norm();311 if (fabs(radius) > MYEPSILON)312 angle = ProjectedVector.Angle(&AngleReferenceVector);313 else314 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues315 316 //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;317 if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)318 {319 angle = 2. * M_PI - angle;320 }321 //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";322 //ProjectedVector.Output(out);323 //*out << endl;324 BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,325 DistancePair (radius, Walker)));326 if (BoundaryTestPair.second)327 { // successfully inserted328 }329 else330 { // same point exists, check first r, then distance of original vectors to center of gravity331 *out << Verbose(2)332 << "Encountered two vectors whose projection onto axis "333 << axis << " is equal: " << endl;334 *out << Verbose(2) << "Present vector: ";335 BoundaryTestPair.first->second.second->x.Output(out);336 *out << endl;337 *out << Verbose(2) << "New vector: ";338 Walker->x.Output(out);339 *out << endl;340 double tmp = ProjectedVector.Norm();341 if (tmp > BoundaryTestPair.first->second.first)342 {343 BoundaryTestPair.first->second.first = tmp;344 BoundaryTestPair.first->second.second = Walker;345 *out << Verbose(2) << "Keeping new vector." << endl;346 }347 else if (tmp == BoundaryTestPair.first->second.first)348 {349 if (BoundaryTestPair.first->second.second->x.ScalarProduct(350 &BoundaryTestPair.first->second.second->x)351 < Walker->x.ScalarProduct(&Walker->x))352 { // Norm() does a sqrt, which makes it a lot slower353 BoundaryTestPair.first->second.second = Walker;354 *out << Verbose(2) << "Keeping new vector." << endl;355 }356 else357 {358 *out << Verbose(2) << "Keeping present vector." << endl;359 }360 }361 else362 {363 *out << Verbose(2) << "Keeping present vector." << endl;364 }365 }366 }367 // printing all inserted for debugging368 //{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 // 3c. throw out points whose distance is less than the mean of left and right neighbours381 bool flag = false;382 do383 { // do as long as we still throw one out per round384 *out << Verbose(1)385 << "Looking for candidates to kick out by convex condition ... "386 << endl;387 flag = false;388 Boundaries::iterator left = BoundaryPoints[axis].end();389 Boundaries::iterator right = BoundaryPoints[axis].end();390 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner391 != BoundaryPoints[axis].end(); runner++)392 {393 // set neighbours correctly394 if (runner == BoundaryPoints[axis].begin())395 {396 left = BoundaryPoints[axis].end();397 }398 else399 {400 left = runner;401 }402 left--;403 right = runner;404 right++;405 if (right == BoundaryPoints[axis].end())406 {407 right = BoundaryPoints[axis].begin();408 }409 // check distance410 411 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)412 {413 Vector SideA, SideB, SideC, SideH;414 SideA.CopyVector(&left->second.second->x);415 SideA.ProjectOntoPlane(&AxisVector);416 //*out << "SideA: ";417 //SideA.Output(out);418 //*out << endl;419 420 SideB.CopyVector(&right->second.second->x);421 SideB.ProjectOntoPlane(&AxisVector);422 //*out << "SideB: ";423 //SideB.Output(out);424 //*out << endl;425 426 SideC.CopyVector(&left->second.second->x);427 SideC.SubtractVector(&right->second.second->x);428 SideC.ProjectOntoPlane(&AxisVector);429 //*out << "SideC: ";430 //SideC.Output(out);431 //*out << endl;432 433 SideH.CopyVector(&runner->second.second->x);434 SideH.ProjectOntoPlane(&AxisVector);435 //*out << "SideH: ";436 //SideH.Output(out);437 //*out << endl;438 439 // calculate each length440 double a = SideA.Norm();441 //double b = SideB.Norm();442 //double c = SideC.Norm();443 double h = SideH.Norm();444 // calculate the angles445 double alpha = SideA.Angle(&SideH);446 double beta = SideA.Angle(&SideC);447 double gamma = SideB.Angle(&SideH);448 double delta = SideC.Angle(&SideH);449 double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha450 < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);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 //*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;453 if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))454 < MYEPSILON) && (h < MinDistance))455 {456 // throw out point457 //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;458 BoundaryPoints[axis].erase(runner);459 flag = true;460 }461 }462 }463 }464 while (flag);465 }466 return BoundaryPoints;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 bool IsAngstroem)481 { 482 // get points on boundary of NULL was given as parameter483 bool BoundaryFreeFlag = false;484 Boundaries *BoundaryPoints = BoundaryPtr;485 if (BoundaryPoints == NULL)486 {487 BoundaryFreeFlag = true;488 BoundaryPoints = GetBoundaryPoints(out, mol);489 }490 else491 {492 *out << Verbose(1) << "Using given boundary points set." << endl;493 }494 // determine biggest "diameter" of cluster for each axis495 Boundaries::iterator Neighbour, OtherNeighbour;496 double *GreatestDiameter = new double[NDIM];497 for (int i = 0; i < NDIM; i++)498 GreatestDiameter[i] = 0.;499 double OldComponent, tmp, w1, w2;500 Vector DistanceVector, OtherVector;501 int component, Othercomponent;502 for (int axis = 0; axis < NDIM; axis++)503 { // regard each projected plane504 //*out << Verbose(1) << "Current axis is " << axis << "." << endl;505 for (int j = 0; j < 2; j++)506 { // and for both axis on the current plane507 component = (axis + j + 1) % NDIM;508 Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;509 //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;510 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner511 != BoundaryPoints[axis].end(); runner++)512 {513 //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;514 // seek for the neighbours pair where the Othercomponent sign flips515 Neighbour = runner;516 Neighbour++;517 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around518 Neighbour = BoundaryPoints[axis].begin();519 DistanceVector.CopyVector(&runner->second.second->x);520 DistanceVector.SubtractVector(&Neighbour->second.second->x);521 do522 { // seek for neighbour pair where it flips523 OldComponent = DistanceVector.x[Othercomponent];524 Neighbour++;525 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around526 Neighbour = BoundaryPoints[axis].begin();527 DistanceVector.CopyVector(&runner->second.second->x);528 DistanceVector.SubtractVector(&Neighbour->second.second->x);529 //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;530 }531 while ((runner != Neighbour) && (fabs(OldComponent / fabs(532 OldComponent) - DistanceVector.x[Othercomponent] / fabs(533 DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip534 if (runner != Neighbour)535 {536 OtherNeighbour = Neighbour;537 if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around538 OtherNeighbour = BoundaryPoints[axis].end();539 OtherNeighbour--;540 //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;541 // now we have found the pair: Neighbour and OtherNeighbour542 OtherVector.CopyVector(&runner->second.second->x);543 OtherVector.SubtractVector(&OtherNeighbour->second.second->x);544 //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;545 //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;546 // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour547 w1 = fabs(OtherVector.x[Othercomponent]);548 w2 = fabs(DistanceVector.x[Othercomponent]);549 tmp = fabs((w1 * DistanceVector.x[component] + w2550 * OtherVector.x[component]) / (w1 + w2));551 // mark if it has greater diameter552 //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;553 GreatestDiameter[component] = (GreatestDiameter[component]554 > tmp) ? GreatestDiameter[component] : tmp;555 } //else556 //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;557 }558 }559 }560 *out << Verbose(0) << "RESULT: The biggest diameters are "561 << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "562 << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"563 : "atomiclength") << "." << endl;564 565 // free reference lists566 if (BoundaryFreeFlag)567 delete[] (BoundaryPoints);568 569 return GreatestDiameter;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 atom *Walker = mol->start;582 bond *Binder = mol->first;583 int i;584 Vector *center = mol->DetermineCenterOfAll(out);585 if (vrmlfile != NULL) {586 //cout << Verbose(1) << "Writing Raster3D file ... ";587 *vrmlfile << "#VRML V2.0 utf8" << endl;588 *vrmlfile << "#Created by molecuilder" << endl;589 *vrmlfile << "#All atoms as spheres" << endl;590 while (Walker->next != mol->end) {591 Walker = Walker->next;592 *vrmlfile << "Sphere {" << endl << ""; // 2 is sphere type593 for (i=0;i<NDIM;i++)594 *vrmlfile << Walker->x.x[i]+center->x[i] << " ";595 *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour596 }597 598 *vrmlfile << "# All bonds as vertices" << endl;599 while (Binder->next != mol->last) {600 Binder = Binder->next;601 *vrmlfile << "3" << endl << ""; // 2 is round-ended cylinder type602 for (i=0;i<NDIM;i++)603 *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";604 *vrmlfile << "\t0.03\t";605 for (i=0;i<NDIM;i++)606 *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";607 *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour608 }609 610 *vrmlfile << "# All tesselation triangles" << endl;611 for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {612 *vrmlfile << "1" << endl << ""; // 1 is triangle type613 for (i=0;i<3;i++) { // print each node614 for (int j=0;j<NDIM;j++)// and for each node all NDIM coordinates615 *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";616 *vrmlfile << "\t";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 } else {622 cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;623 }624 delete(center);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 atom *Walker = mol->start;636 bond *Binder = mol->first;637 int i;638 Vector *center = mol->DetermineCenterOfAll(out);639 if (rasterfile != NULL) {640 //cout << Verbose(1) << "Writing Raster3D file ... ";641 *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;642 *rasterfile << "@header.r3d" << endl;643 *rasterfile << "# All atoms as spheres" << endl;644 while (Walker->next != mol->end) {645 Walker = Walker->next;646 *rasterfile << "2" << endl << " ";// 2 is sphere type647 for (i=0;i<NDIM;i++)648 *rasterfile << Walker->x.x[i]+center->x[i] << " ";649 *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour650 }651 652 *rasterfile << "# All bonds as vertices" << endl;653 while (Binder->next != mol->last) {654 Binder = Binder->next;655 *rasterfile << "3" << endl << " ";// 2 is round-ended cylinder type656 for (i=0;i<NDIM;i++)657 *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";658 *rasterfile << "\t0.03\t";659 for (i=0;i<NDIM;i++)660 *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";661 *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour662 }663 664 *rasterfile << "# All tesselation triangles" << endl;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 for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {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 *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";671 *rasterfile << "\t";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 } else {678 cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;679 }680 delete(center);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 class Tesselation *TesselStruct, class molecule *mol, int N)691 { 692 if (tecplot != NULL)693 {694 *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;695 *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;696 *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N="697 << TesselStruct->PointsOnBoundaryCount << ", E="698 << TesselStruct->TrianglesOnBoundaryCount699 << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;700 int *LookupList = new int[mol->AtomCount];701 for (int i = 0; i < mol->AtomCount; i++)702 LookupList[i] = -1;703 704 // print atom coordinates705 *out << Verbose(2) << "The following triangles were created:";706 int Counter = 1;707 atom *Walker = NULL;708 for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target709 != TesselStruct->PointsOnBoundary.end(); target++)710 {711 Walker = target->second->node;712 LookupList[Walker->nr] = Counter++;713 *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " "714 << Walker->x.x[2] << " " << endl;715 }716 *tecplot << endl;717 // print connectivity718 for (TriangleMap::iterator runner =719 TesselStruct->TrianglesOnBoundary.begin(); runner720 != TesselStruct->TrianglesOnBoundary.end(); runner++)721 {722 *out << " " << runner->second->endpoints[0]->node->Name << "<->"723 << runner->second->endpoints[1]->node->Name << "<->"724 << runner->second->endpoints[2]->node->Name;725 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " "726 << LookupList[runner->second->endpoints[1]->node->nr] << " "727 << LookupList[runner->second->endpoints[2]->node->nr] << endl;728 }729 delete[] (LookupList);730 *out << endl;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 Boundaries *BoundaryPtr, molecule *mol)746 { 747 bool IsAngstroem = configuration->GetIsAngstroem();748 atom *Walker = NULL;749 struct Tesselation *TesselStruct = new Tesselation;750 bool BoundaryFreeFlag = false;751 Boundaries *BoundaryPoints = BoundaryPtr;752 double volume = 0.;753 double PyramidVolume = 0.;754 double G, h;755 Vector x, y;756 double a, b, c;757 758 //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.759 760 // 1. calculate center of gravity761 *out << endl;762 Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);763 764 // 2. translate all points into CoG765 *out << Verbose(1) << "Translating system to Center of Gravity." << endl;766 Walker = mol->start;767 while (Walker->next != mol->end)768 {769 Walker = Walker->next;770 Walker->x.Translate(CenterOfGravity);771 }772 773 // 3. Find all points on the boundary774 if (BoundaryPoints == NULL)775 {776 BoundaryFreeFlag = true;777 BoundaryPoints = GetBoundaryPoints(out, mol);778 }779 else780 {781 *out << Verbose(1) << "Using given boundary points set." << endl;782 }783 784 // 4. fill the boundary point list785 for (int axis = 0; axis < NDIM; axis++)786 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner787 != BoundaryPoints[axis].end(); runner++)788 {789 TesselStruct->AddPoint(runner->second.second);790 }791 792 *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount793 << " points on the convex boundary." << endl;794 // now we have the whole set of edge points in the BoundaryList795 796 // listing for debugging797 //*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 // 5a. guess starting triangle804 TesselStruct->GuessStartingTriangle(out);805 806 // 5b. go through all lines, that are not yet part of two triangles (only of one so far)807 TesselStruct->TesselateOnBoundary(out, configuration, mol);808 809 *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount810 << " triangles with " << TesselStruct->LinesOnBoundaryCount811 << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."812 << endl;813 814 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes815 *out << Verbose(1)816 << "Calculating the volume of the pyramids formed out of triangles and center of gravity."817 << endl;818 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner819 != TesselStruct->TrianglesOnBoundary.end(); runner++)820 { // go through every triangle, calculate volume of its pyramid with CoG as peak821 x.CopyVector(&runner->second->endpoints[0]->node->x);822 x.SubtractVector(&runner->second->endpoints[1]->node->x);823 y.CopyVector(&runner->second->endpoints[0]->node->x);824 y.SubtractVector(&runner->second->endpoints[2]->node->x);825 a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(826 &runner->second->endpoints[1]->node->x));827 b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(828 &runner->second->endpoints[2]->node->x));829 c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(830 &runner->second->endpoints[1]->node->x));831 G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle832 x.MakeNormalVector(&runner->second->endpoints[0]->node->x,833 &runner->second->endpoints[1]->node->x,834 &runner->second->endpoints[2]->node->x);835 x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));836 h = x.Norm(); // distance of CoG to triangle837 PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)838 *out << Verbose(2) << "Area of triangle is " << G << " "839 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "840 << h << " and the volume is " << PyramidVolume << " "841 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;842 volume += PyramidVolume;843 }844 *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)845 << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."846 << endl;847 848 // 7. translate all points back from CoG849 *out << Verbose(1) << "Translating system back from Center of Gravity."850 << endl;851 CenterOfGravity->Scale(-1);852 Walker = mol->start;853 while (Walker->next != mol->end)854 {855 Walker = Walker->next;856 Walker->x.Translate(CenterOfGravity);857 }858 859 // 8. Store triangles in tecplot file860 string OutputName(filename);861 OutputName.append(TecplotSuffix);862 ofstream *tecplot = new ofstream(OutputName.c_str());863 write_tecplot_file(out, tecplot, TesselStruct, mol, 0);864 tecplot->close();865 delete(tecplot);866 867 // free reference lists868 if (BoundaryFreeFlag)869 delete[] (BoundaryPoints);870 871 return volume;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 double ClusterVolume, double celldensity)886 { 887 // transform to PAS888 mol->PrincipalAxisSystem(out, true);889 890 // some preparations beforehand891 bool IsAngstroem = configuration->GetIsAngstroem();892 Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);893 double clustervolume;894 if (ClusterVolume == 0)895 clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,896 BoundaryPoints, mol);897 else898 clustervolume = ClusterVolume;899 double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,900 IsAngstroem);901 Vector BoxLengths;902 int repetition[NDIM] =903 { 1, 1, 1 };904 int TotalNoClusters = 1;905 for (int i = 0; i < NDIM; i++)906 TotalNoClusters *= repetition[i];907 908 // sum up the atomic masses909 double totalmass = 0.;910 atom *Walker = mol->start;911 while (Walker->next != mol->end)912 {913 Walker = Walker->next;914 totalmass += Walker->type->mass;915 }916 *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)917 << totalmass << " atomicmassunit." << endl;918 919 *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)920 << totalmass / clustervolume << " atomicmassunit/"921 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;922 923 // solve cubic polynomial924 *out << Verbose(1) << "Solving equidistant suspension in water problem ..."925 << endl;926 double cellvolume;927 if (IsAngstroem)928 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass929 / clustervolume)) / (celldensity - 1);930 else931 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass932 / clustervolume)) / (celldensity - 1);933 *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity934 << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"935 : "atomiclength") << "^3." << endl;936 937 double minimumvolume = TotalNoClusters * (GreatestDiameter[0]938 * GreatestDiameter[1] * GreatestDiameter[2]);939 *out << Verbose(1)940 << "Minimum volume of the convex envelope contained in a rectangular box is "941 << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"942 : "atomiclength") << "^3." << endl;943 if (minimumvolume > cellvolume)944 {945 cerr << Verbose(0)946 << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"947 << endl;948 cout << Verbose(0)949 << "Setting Box dimensions to minimum possible, the greatest diameters."950 << endl;951 for (int i = 0; i < NDIM; i++)952 BoxLengths.x[i] = GreatestDiameter[i];953 mol->CenterEdge(out, &BoxLengths);954 }955 else956 {957 BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]958 * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);959 BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]960 * GreatestDiameter[1] + repetition[0] * repetition[2]961 * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]962 * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);963 BoxLengths.x[2] = minimumvolume - cellvolume;964 double x0 = 0., x1 = 0., x2 = 0.;965 if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],966 BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return967 *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0968 << " ." << endl;969 else970 {971 *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0972 << " and " << x1 << " and " << x2 << " ." << endl;973 x0 = x2; // sorted in ascending order974 }975 976 cellvolume = 1;977 for (int i = 0; i < NDIM; i++)978 {979 BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);980 cellvolume *= BoxLengths.x[i];981 }982 983 // set new box dimensions984 *out << Verbose(0) << "Translating to box with these boundaries." << endl;985 mol->CenterInBox((ofstream *) &cout, &BoxLengths);986 }987 // update Box of atoms by boundary988 mol->SetBoxDimension(&BoxLengths);989 *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "990 << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "991 << BoxLengths.x[2] << " with total volume of " << cellvolume << " "992 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;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 PointsOnBoundaryCount = 0;1003 LinesOnBoundaryCount = 0;1004 TrianglesOnBoundaryCount = 0;1005 TriangleFilesWritten = 0;1023 PointsOnBoundaryCount = 0; 1024 LinesOnBoundaryCount = 0; 1025 TrianglesOnBoundaryCount = 0; 1026 TriangleFilesWritten = 0; 1006 1027 } 1007 1028 ; … … 1012 1033 Tesselation::~Tesselation() 1013 1034 { 1014 cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;1015 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {1016 if (runner->second != NULL) {1017 delete (runner->second);1018 runner->second = NULL;1019 } else1020 cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;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 // 4b. create a starting triangle1034 // 4b1. create all distances1035 DistanceMultiMap DistanceMMap;1036 double distance, tmp;1037 Vector PlaneVector, TrialVector;1038 PointMap::iterator A, B, C; // three nodes of the first triangle1039 A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily1040 1041 // with A chosen, take each pair B,C and sort1042 if (A != PointsOnBoundary.end())1043 {1044 B = A;1045 B++;1046 for (; B != PointsOnBoundary.end(); B++)1047 {1048 C = B;1049 C++;1050 for (; C != PointsOnBoundary.end(); C++)1051 {1052 tmp = A->second->node->x.DistanceSquared(&B->second->node->x);1053 distance = tmp * tmp;1054 tmp = A->second->node->x.DistanceSquared(&C->second->node->x);1055 distance += tmp * tmp;1056 tmp = B->second->node->x.DistanceSquared(&C->second->node->x);1057 distance += tmp * tmp;1058 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<1059 PointMap::iterator, PointMap::iterator> (B, C)));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 // 4b2. pick three baselines forming a triangle1070 // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate1071 DistanceMultiMap::iterator baseline = DistanceMMap.begin();1072 for (; baseline != DistanceMMap.end(); baseline++)1073 {1074 // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate1075 // 2. next, we have to check whether all points reside on only one side of the triangle1076 // 3. construct plane vector1077 PlaneVector.MakeNormalVector(&A->second->node->x,1078 &baseline->second.first->second->node->x,1079 &baseline->second.second->second->node->x);1080 *out << Verbose(2) << "Plane vector of candidate triangle is ";1081 PlaneVector.Output(out);1082 *out << endl;1083 // 4. loop over all points1084 double sign = 0.;1085 PointMap::iterator checker = PointsOnBoundary.begin();1086 for (; checker != PointsOnBoundary.end(); checker++)1087 {1088 // (neglecting A,B,C)1089 if ((checker == A) || (checker == baseline->second.first) || (checker1090 == baseline->second.second))1091 continue;1092 // 4a. project onto plane vector1093 TrialVector.CopyVector(&checker->second->node->x);1094 TrialVector.SubtractVector(&A->second->node->x);1095 distance = TrialVector.Projection(&PlaneVector);1096 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok1097 continue;1098 *out << Verbose(3) << "Projection of " << checker->second->node->Name1099 << " yields distance of " << distance << "." << endl;1100 tmp = distance / fabs(distance);1101 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)1102 if ((sign != 0) && (tmp != sign))1103 {1104 // 4c. If so, break 4. loop and continue with next candidate in 1. loop1105 *out << Verbose(2) << "Current candidates: "1106 << A->second->node->Name << ","1107 << baseline->second.first->second->node->Name << ","1108 << baseline->second.second->second->node->Name << " leave "1109 << checker->second->node->Name << " outside the convex hull."1110 << endl;1111 break;1112 }1113 else1114 { // note the sign for later1115 *out << Verbose(2) << "Current candidates: "1116 << A->second->node->Name << ","1117 << baseline->second.first->second->node->Name << ","1118 << baseline->second.second->second->node->Name << " leave "1119 << checker->second->node->Name << " inside the convex hull."1120 << endl;1121 sign = tmp;1122 }1123 // 4d. Check whether the point is inside the triangle (check distance to each node1124 tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);1125 int innerpoint = 0;1126 if ((tmp < A->second->node->x.DistanceSquared(1127 &baseline->second.first->second->node->x)) && (tmp1128 < A->second->node->x.DistanceSquared(1129 &baseline->second.second->second->node->x)))1130 innerpoint++;1131 tmp = checker->second->node->x.DistanceSquared(1132 &baseline->second.first->second->node->x);1133 if ((tmp < baseline->second.first->second->node->x.DistanceSquared(1134 &A->second->node->x)) && (tmp1135 < baseline->second.first->second->node->x.DistanceSquared(1136 &baseline->second.second->second->node->x)))1137 innerpoint++;1138 tmp = checker->second->node->x.DistanceSquared(1139 &baseline->second.second->second->node->x);1140 if ((tmp < baseline->second.second->second->node->x.DistanceSquared(1141 &baseline->second.first->second->node->x)) && (tmp1142 < baseline->second.second->second->node->x.DistanceSquared(1143 &A->second->node->x)))1144 innerpoint++;1145 // 4e. If so, break 4. loop and continue with next candidate in 1. loop1146 if (innerpoint == 3)1147 break;1148 }1149 // 5. come this far, all on same side? Then break 1. loop and construct triangle1150 if (checker == PointsOnBoundary.end())1151 {1152 *out << "Looks like we have a candidate!" << endl;1153 break;1154 }1155 }1156 if (baseline != DistanceMMap.end())1157 {1158 BPS[0] = baseline->second.first->second;1159 BPS[1] = baseline->second.second->second;1160 BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);1161 BPS[0] = A->second;1162 BPS[1] = baseline->second.second->second;1163 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);1164 BPS[0] = baseline->second.first->second;1165 BPS[1] = A->second;1166 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);1167 1168 // 4b3. insert created triangle1169 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);1170 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));1171 TrianglesOnBoundaryCount++;1172 for (int i = 0; i < NDIM; i++)1173 {1174 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));1175 LinesOnBoundaryCount++;1176 }1177 1178 *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;1179 }1180 else1181 {1182 *out << Verbose(1) << "No starting triangle found." << endl;1183 exit(255);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 * -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors1194 * -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to1195 * baseline in triangle plane pointing out of the triangle and normal vector of new triangle)1196 * -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)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 molecule *mol)1204 { 1205 bool flag;1206 PointMap::iterator winner;1207 class BoundaryPointSet *peak = NULL;1208 double SmallestAngle, TempAngle;1209 Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,1210 PropagationVector;1211 LineMap::iterator LineChecker[2];1212 do1213 {1214 flag = false;1215 for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline1216 != LinesOnBoundary.end(); baseline++)1217 if (baseline->second->TrianglesCount == 1)1218 {1219 *out << Verbose(2) << "Current baseline is between "1220 << *(baseline->second) << "." << endl;1221 // 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)1222 SmallestAngle = M_PI;1223 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far1224 // get peak point with respect to this base line's only triangle1225 for (int i = 0; i < 3; i++)1226 if ((BTS->endpoints[i] != baseline->second->endpoints[0])1227 && (BTS->endpoints[i] != baseline->second->endpoints[1]))1228 peak = BTS->endpoints[i];1229 *out << Verbose(3) << " and has peak " << *peak << "." << endl;1230 // normal vector of triangle1231 BTS->GetNormalVector(NormalVector);1232 *out << Verbose(4) << "NormalVector of base triangle is ";1233 NormalVector.Output(out);1234 *out << endl;1235 // offset to center of triangle1236 CenterVector.Zero();1237 for (int i = 0; i < 3; i++)1238 CenterVector.AddVector(&BTS->endpoints[i]->node->x);1239 CenterVector.Scale(1. / 3.);1240 *out << Verbose(4) << "CenterVector of base triangle is ";1241 CenterVector.Output(out);1242 *out << endl;1243 // vector in propagation direction (out of triangle)1244 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)1245 TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);1246 TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);1247 PropagationVector.MakeNormalVector(&TempVector, &NormalVector);1248 TempVector.CopyVector(&CenterVector);1249 TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!1250 //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;1251 if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline1252 PropagationVector.Scale(-1.);1253 *out << Verbose(4) << "PropagationVector of base triangle is ";1254 PropagationVector.Output(out);1255 *out << endl;1256 winner = PointsOnBoundary.end();1257 for (PointMap::iterator target = PointsOnBoundary.begin(); target1258 != PointsOnBoundary.end(); target++)1259 if ((target->second != baseline->second->endpoints[0])1260 && (target->second != baseline->second->endpoints[1]))1261 { // don't take the same endpoints1262 *out << Verbose(3) << "Target point is " << *(target->second)1263 << ":";1264 bool continueflag = true;1265 1266 VirtualNormalVector.CopyVector(1267 &baseline->second->endpoints[0]->node->x);1268 VirtualNormalVector.AddVector(1269 &baseline->second->endpoints[0]->node->x);1270 VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line1271 VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target1272 TempAngle = VirtualNormalVector.Angle(&PropagationVector);1273 continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)1274 if (!continueflag)1275 {1276 *out << Verbose(4)1277 << "Angle between propagation direction and base line to "1278 << *(target->second) << " is " << TempAngle1279 << ", bad direction!" << endl;1280 continue;1281 }1282 else1283 *out << Verbose(4)1284 << "Angle between propagation direction and base line to "1285 << *(target->second) << " is " << TempAngle1286 << ", good direction!" << endl;1287 LineChecker[0] = baseline->second->endpoints[0]->lines.find(1288 target->first);1289 LineChecker[1] = baseline->second->endpoints[1]->lines.find(1290 target->first);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 // check first endpoint (if any connecting line goes to target or at least not more than 1)1300 continueflag = continueflag && (((LineChecker[0]1301 == baseline->second->endpoints[0]->lines.end())1302 || (LineChecker[0]->second->TrianglesCount == 1)));1303 if (!continueflag)1304 {1305 *out << Verbose(4) << *(baseline->second->endpoints[0])1306 << " has line " << *(LineChecker[0]->second)1307 << " to " << *(target->second)1308 << " as endpoint with "1309 << LineChecker[0]->second->TrianglesCount1310 << " triangles." << endl;1311 continue;1312 }1313 // check second endpoint (if any connecting line goes to target or at least not more than 1)1314 continueflag = continueflag && (((LineChecker[1]1315 == baseline->second->endpoints[1]->lines.end())1316 || (LineChecker[1]->second->TrianglesCount == 1)));1317 if (!continueflag)1318 {1319 *out << Verbose(4) << *(baseline->second->endpoints[1])1320 << " has line " << *(LineChecker[1]->second)1321 << " to " << *(target->second)1322 << " as endpoint with "1323 << LineChecker[1]->second->TrianglesCount1324 << " triangles." << endl;1325 continue;1326 }1327 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)1328 continueflag = continueflag && (!(((LineChecker[0]1329 != baseline->second->endpoints[0]->lines.end())1330 && (LineChecker[1]1331 != baseline->second->endpoints[1]->lines.end())1332 && (GetCommonEndpoint(LineChecker[0]->second,1333 LineChecker[1]->second) == peak))));1334 if (!continueflag)1335 {1336 *out << Verbose(4) << "Current target is peak!" << endl;1337 continue;1338 }1339 // in case NOT both were found1340 if (continueflag)1341 { // create virtually this triangle, get its normal vector, calculate angle1342 flag = true;1343 VirtualNormalVector.MakeNormalVector(1344 &baseline->second->endpoints[0]->node->x,1345 &baseline->second->endpoints[1]->node->x,1346 &target->second->node->x);1347 // make it always point inward1348 if (baseline->second->endpoints[0]->node->x.Projection(1349 &VirtualNormalVector) > 0)1350 VirtualNormalVector.Scale(-1.);1351 // calculate angle1352 TempAngle = NormalVector.Angle(&VirtualNormalVector);1353 *out << Verbose(4) << "NormalVector is ";1354 VirtualNormalVector.Output(out);1355 *out << " and the angle is " << TempAngle << "." << endl;1356 if (SmallestAngle > TempAngle)1357 { // set to new possible winner1358 SmallestAngle = TempAngle;1359 winner = target;1360 }1361 }1362 }1363 // 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 triangle1364 if (winner != PointsOnBoundary.end())1365 {1366 *out << Verbose(2) << "Winning target point is "1367 << *(winner->second) << " with angle " << SmallestAngle1368 << "." << endl;1369 // create the lins of not yet present1370 BLS[0] = baseline->second;1371 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)1372 LineChecker[0] = baseline->second->endpoints[0]->lines.find(1373 winner->first);1374 LineChecker[1] = baseline->second->endpoints[1]->lines.find(1375 winner->first);1376 if (LineChecker[0]1377 == baseline->second->endpoints[0]->lines.end())1378 { // create1379 BPS[0] = baseline->second->endpoints[0];1380 BPS[1] = winner->second;1381 BLS[1] = new class BoundaryLineSet(BPS,1382 LinesOnBoundaryCount);1383 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,1384 BLS[1]));1385 LinesOnBoundaryCount++;1386 }1387 else1388 BLS[1] = LineChecker[0]->second;1389 if (LineChecker[1]1390 == baseline->second->endpoints[1]->lines.end())1391 { // create1392 BPS[0] = baseline->second->endpoints[1];1393 BPS[1] = winner->second;1394 BLS[2] = new class BoundaryLineSet(BPS,1395 LinesOnBoundaryCount);1396 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,1397 BLS[2]));1398 LinesOnBoundaryCount++;1399 }1400 else1401 BLS[2] = LineChecker[1]->second;1402 BTS = new class BoundaryTriangleSet(BLS,1403 TrianglesOnBoundaryCount);1404 TrianglesOnBoundary.insert(TrianglePair(1405 TrianglesOnBoundaryCount, BTS));1406 TrianglesOnBoundaryCount++;1407 }1408 else1409 {1410 *out << Verbose(1)1411 << "I could not determine a winner for this baseline "1412 << *(baseline->second) << "." << endl;1413 }1414 1415 // 5d. If the set of lines is not yet empty, go to 5. and continue1416 }1417 else1418 *out << Verbose(2) << "Baseline candidate " << *(baseline->second)1419 << " has a triangle count of "1420 << baseline->second->TrianglesCount << "." << endl;1421 }1422 while (flag);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 PointTestPair InsertUnique;1434 BPS[0] = new class BoundaryPointSet(Walker);1435 InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));1436 if (InsertUnique.second) // if new point was not present before, increase counter1437 PointsOnBoundaryCount++;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 PointTestPair InsertUnique;1450 TPS[n] = new class BoundaryPointSet(Candidate);1451 InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));1452 if (InsertUnique.second) // if new point was not present before, increase counter1453 {1454 PointsOnBoundaryCount++;1455 }1456 else1457 {1458 delete TPS[n];1459 cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)1460 << " gibt's schon in der PointMap." << endl;1461 TPS[n] = (InsertUnique.first)->second;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 cout << Verbose(1) << "Adding triangle to its lines" << endl;1523 int i = 0;1524 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));1525 TrianglesOnBoundaryCount++;1526 1527 /*1528 * this is apparently done when constructing triangle1529 1530 for (i=0; i<3; i++)1531 {1532 BLS[i]->AddTriangle(BTS);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 inPlace = 1 => A is replaced with the LU decomposed copy.1542 inPlace = 0 => A is retained, and a copy is used for LU.1543 */1544 1545 double det;1546 int signum;1547 gsl_permutation *p = gsl_permutation_alloc(A->size1);1548 gsl_matrix *tmpA;1549 1550 if (inPlace)1551 tmpA = A;1552 else {1553 gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);1554 gsl_matrix_memcpy(tmpA , A);1555 }1556 1557 1558 gsl_linalg_LU_decomp(tmpA , p , &signum);1559 det = gsl_linalg_LU_det(tmpA , signum);1560 gsl_permutation_free(p);1561 if (! inPlace)1562 gsl_matrix_free(tmpA);1563 1564 return det;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 gsl_matrix *A = gsl_matrix_calloc(3,3);1570 double m11, m12, m13, m14;1571 1572 for(int i=0;i<3;i++) {1573 gsl_matrix_set(A, i, 0, a.x[i]);1574 gsl_matrix_set(A, i, 1, b.x[i]);1575 gsl_matrix_set(A, i, 2, c.x[i]);1576 }1577 m11 = det_get(A, 1);1578 1579 for(int i=0;i<3;i++) {1580 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);1581 gsl_matrix_set(A, i, 1, b.x[i]);1582 gsl_matrix_set(A, i, 2, c.x[i]);1583 }1584 m12 = det_get(A, 1);1585 1586 for(int i=0;i<3;i++) {1587 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);1588 gsl_matrix_set(A, i, 1, a.x[i]);1589 gsl_matrix_set(A, i, 2, c.x[i]);1590 }1591 m13 = det_get(A, 1);1592 1593 for(int i=0;i<3;i++) {1594 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);1595 gsl_matrix_set(A, i, 1, a.x[i]);1596 gsl_matrix_set(A, i, 2, b.x[i]);1597 }1598 m14 = det_get(A, 1);1599 1600 if (fabs(m11) < MYEPSILON)1601 cerr << "ERROR: three points are colinear." << endl;1602 1603 center->x[0] =0.5 * m12/ m11;1604 center->x[1] = -0.5 * m13/ m11;1605 center->x[2] =0.5 * m14/ m11;1606 1607 if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)1608 cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;1609 1610 gsl_matrix_free(A);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 double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)1633 { 1634 Vector TempNormal, helper;1635 double Restradius;1636 Vector OtherCenter;1637 cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";1638 Center->Zero();1639 helper.CopyVector(&a);1640 helper.Scale(sin(2.*alpha));1641 Center->AddVector(&helper);1642 helper.CopyVector(&b);1643 helper.Scale(sin(2.*beta));1644 Center->AddVector(&helper);1645 helper.CopyVector(&c);1646 helper.Scale(sin(2.*gamma));1647 Center->AddVector(&helper);1648 //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;1649 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));1650 NewUmkreismittelpunkt->CopyVector(Center);1651 cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";1652 // Here we calculated center of circumscribing circle, using barycentric coordinates1653 cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";1654 1655 TempNormal.CopyVector(&a);1656 TempNormal.SubtractVector(&b);1657 helper.CopyVector(&a);1658 helper.SubtractVector(&c);1659 TempNormal.VectorProduct(&helper);1660 if (fabs(HalfplaneIndicator) < MYEPSILON)1661 {1662 if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))1663 {1664 TempNormal.Scale(-1);1665 }1666 }1667 else1668 {1669 if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)1670 {1671 TempNormal.Scale(-1);1672 }1673 }1674 1675 TempNormal.Normalize();1676 Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);1677 cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";1678 TempNormal.Scale(Restradius);1679 cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";1680 1681 Center->AddVector(&TempNormal);1682 cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";1683 get_sphere(&OtherCenter, a, b, c, RADIUS);1684 cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";1685 cout << Verbose(3) << "End of Get_center_of_sphere.\n";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 * supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \1691 * upon which we operate.1692 * If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \1693 * direction and angle into Storage.1694 * We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \1695 * with all neighbours of the candidate.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 Vector helper;1877 double alpha, beta, gamma;1878 Vector SideA, SideB, SideC;1879 SideA.CopyVector(b);1880 SideA.SubtractVector(c);1881 SideB.CopyVector(c);1882 SideB.SubtractVector(a);1883 SideC.CopyVector(a);1884 SideC.SubtractVector(b);1885 alpha = M_PI - SideB.Angle(&SideC);1886 beta = M_PI - SideC.Angle(&SideA);1887 gamma = M_PI - SideA.Angle(&SideB);1888 cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;1889 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)1890 cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;1891 1892 Center->Zero();1893 helper.CopyVector(a);1894 helper.Scale(sin(2.*alpha));1895 Center->AddVector(&helper);1896 helper.CopyVector(b);1897 helper.Scale(sin(2.*beta));1898 Center->AddVector(&helper);1899 helper.CopyVector(c);1900 helper.Scale(sin(2.*gamma));1901 Center->AddVector(&helper);1902 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));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 Vector helper;1921 double radius, alpha;1922 1923 helper.CopyVector(&NewSphereCenter);1924 // test whether new center is on the parameter circle's plane1925 if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {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 helper.ProjectOntoPlane(&CirclePlaneNormal);1928 }1929 radius = helper.ScalarProduct(&helper);1930 // test whether the new center vector has length of CircleRadius1931 if (fabs(radius - CircleRadius) > HULLEPSILON)1932 cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;1933 alpha = helper.Angle(&OldSphereCenter);1934 // make the angle unique by checking the halfplanes/search direction1935 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals1936 alpha = 2.*M_PI - alpha;1937 cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;1938 radius = helper.Distance(&OldSphereCenter);1939 helper.ProjectOntoPlane(&NormalVector);1940 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles1941 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {1942 cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;1943 return alpha;1944 } else {1945 cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;1946 return 2.*M_PI;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 // atom *Walker = NULL;1983 // Vector CircleCenter;// center of the circle, i.e. of the band of sphere's centers1984 // Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in1985 // 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 // Vector helper;1991 // LinkedAtoms *List = NULL;1992 // double CircleRadius; // radius of this circle1993 // double radius;1994 // double alpha, Otheralpha; // angles (i.e. parameter for the circle).1995 // double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle1996 // int N[NDIM], Nlower[NDIM], Nupper[NDIM];1997 // atom *Candidate = NULL;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 // cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;2038 // cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl; 2000 2039 // 2001 // cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;2040 // cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl; 2002 2041 // 2003 // // construct center of circle2004 // CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));2005 // CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);2006 // CircleCenter.Scale(0.5);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 // // construct normal vector of circle2009 // CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);2010 // CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);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 // // calculate squared radius of circle2013 // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);2014 // if (radius/4. < RADIUS*RADIUS) {2015 // CircleRadius = RADIUS*RADIUS - radius/4.;2016 // CirclePlaneNormal.Normalize();2017 // cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;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 // // construct old center2020 // GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));2021 // helper.CopyVector(&BaseTriangle->NormalVector);// normal vector ensures that this is correct center of the two possible ones2022 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);2023 // helper.Scale(sqrt(RADIUS*RADIUS - radius));2024 // OldSphereCenter.AddVector(&helper);2025 // OldSphereCenter.SubtractVector(&CircleCenter);2026 // cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;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 // // test whether old center is on the band's plane2029 // if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {2030 // cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;2031 // OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);2032 // }2033 // radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);2034 // if (fabs(radius - CircleRadius) < HULLEPSILON) {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 // // construct SearchDirection2037 // SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);2038 // helper.CopyVector(&BaseLine->endpoints[0]->node->x);2039 // for(int i=0;i<3;i++)// just take next different endpoint2040 // if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {2041 // helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);2042 // }2043 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!2044 // SearchDirection.Scale(-1.);2045 // SearchDirection.ProjectOntoPlane(&OldSphereCenter);2046 // SearchDirection.Normalize();2047 // cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;2048 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {// rotated the wrong way!2049 // cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;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 // for(int i=0;i<NDIM;i++) // store indices of this cell2054 // N[i] = LC->n[i];2055 // cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;2056 // } else {2057 // cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;2058 // return;2059 // }2060 // // then go through the current and all neighbouring cells and check the contained atoms for possible candidates2061 // cout << Verbose(2) << "LC Intervals:";2062 // for (int i=0;i<NDIM;i++) {2063 // Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;2064 // Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;2065 // cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";2066 // }2067 // cout << endl;2068 // for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)2069 // for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)2070 // for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {2071 // List = LC->GetCurrentCell();2072 // cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;2073 // if (List != NULL) {2074 // for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {2075 // Candidate = (*Runner);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 // // check for three unique points2078 // if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {2079 // cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;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 // // construct both new centers2082 // GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));2083 // OtherNewSphereCenter.CopyVector(&NewSphereCenter);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 // if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {2086 // helper.CopyVector(&NewNormalVector);2087 // cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;2088 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);2089 // if (radius < RADIUS*RADIUS) {2090 // helper.Scale(sqrt(RADIUS*RADIUS - radius));2091 // cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;2092 // NewSphereCenter.AddVector(&helper);2093 // NewSphereCenter.SubtractVector(&CircleCenter);2094 // cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;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 // helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction2097 // OtherNewSphereCenter.AddVector(&helper);2098 // OtherNewSphereCenter.SubtractVector(&CircleCenter);2099 // cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;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 // // check both possible centers2102 // alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);2103 // Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);2104 // alpha = min(alpha, Otheralpha);2105 // if (*ShortestAngle > alpha) {2106 // OptCandidate = Candidate;2107 // *ShortestAngle = alpha;2108 // if (alpha != Otheralpha)2109 // OptCandidateCenter->CopyVector(&NewSphereCenter);2110 // else2111 // OptCandidateCenter->CopyVector(&OtherNewSphereCenter);2112 // cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;2113 // } else {2114 // if (OptCandidate != NULL)2115 // cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;2116 // else2117 // cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;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 // } else {2121 // cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;2122 // }2123 // } else {2124 // cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;2125 // }2126 // } else {2127 // cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;2128 // }2129 // }2130 // }2131 // }2132 // } else {2133 // cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;2134 // }2135 // } else {2136 // cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;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 // cout << Verbose(1) << "End of Find_next_suitable_point" << endl;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 if (endpoint1->Nr < endpoint2->Nr) {30 endpoints[0] = endpoint1;31 endpoints[1] = endpoint2;32 } else {33 endpoints[0] = endpoint2;34 endpoints[1] = endpoint1;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 public:40 BoundaryPointSet();41 BoundaryPointSet(atom *Walker);42 ~BoundaryPointSet();36 public: 37 BoundaryPointSet(); 38 BoundaryPointSet(atom *Walker); 39 ~BoundaryPointSet(); 43 40 44 void AddLine(class BoundaryLineSet *line);41 void AddLine(class BoundaryLineSet *line); 45 42 46 LineMap lines;47 int LinesCount;48 atom *node;49 int Nr;43 LineMap lines; 44 int LinesCount; 45 atom *node; 46 int Nr; 50 47 }; 51 48 52 49 class BoundaryLineSet { 53 public:54 BoundaryLineSet();55 BoundaryLineSet(class BoundaryPointSet *Point[2], int number);56 ~BoundaryLineSet();50 public: 51 BoundaryLineSet(); 52 BoundaryLineSet(class BoundaryPointSet *Point[2], int number); 53 ~BoundaryLineSet(); 57 54 58 void AddTriangle(class BoundaryTriangleSet *triangle);55 void AddTriangle(class BoundaryTriangleSet *triangle); 59 56 60 class BoundaryPointSet *endpoints[2];61 TriangleMap triangles;62 int TrianglesCount;63 int Nr;57 class BoundaryPointSet *endpoints[2]; 58 TriangleMap triangles; 59 int TrianglesCount; 60 int Nr; 64 61 }; 65 62 66 63 class BoundaryTriangleSet { 67 public:68 BoundaryTriangleSet();69 BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);70 ~BoundaryTriangleSet();64 public: 65 BoundaryTriangleSet(); 66 BoundaryTriangleSet(class BoundaryLineSet *line[3], int number); 67 ~BoundaryTriangleSet(); 71 68 72 void GetNormalVector(Vector &NormalVector);69 void GetNormalVector(Vector &NormalVector); 73 70 74 class BoundaryPointSet *endpoints[3];75 class BoundaryLineSet *lines[3];76 Vector NormalVector;77 int Nr;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 public:90 public: 82 91 83 Tesselation();84 ~Tesselation();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 PointMap PointsOnBoundary;98 LineMap LinesOnBoundary;99 TriangleMap TrianglesOnBoundary;100 class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions101 class BoundaryPointSet *BPS[2];102 class BoundaryLineSet *BLS[3];103 class BoundaryTriangleSet *BTS;104 int PointsOnBoundaryCount;105 int LinesOnBoundaryCount;106 int TrianglesOnBoundaryCount;107 int TriangleFilesWritten;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.
