| 1 | #include "boundary.hpp"
|
|---|
| 2 |
|
|---|
| 3 | #define DEBUG 1
|
|---|
| 4 | #define DoSingleStepOutput 0
|
|---|
| 5 | #define DoTecplotOutput 1
|
|---|
| 6 | #define DoRaster3DOutput 1
|
|---|
| 7 | #define DoVRMLOutput 1
|
|---|
| 8 | #define TecplotSuffix ".dat"
|
|---|
| 9 | #define Raster3DSuffix ".r3d"
|
|---|
| 10 | #define VRMLSUffix ".wrl"
|
|---|
| 11 | #define HULLEPSILON MYEPSILON
|
|---|
| 12 |
|
|---|
| 13 | // ======================================== Points on Boundary =================================
|
|---|
| 14 |
|
|---|
| 15 | BoundaryPointSet::BoundaryPointSet()
|
|---|
| 16 | {
|
|---|
| 17 | LinesCount = 0;
|
|---|
| 18 | Nr = -1;
|
|---|
| 19 | }
|
|---|
| 20 | ;
|
|---|
| 21 |
|
|---|
| 22 | BoundaryPointSet::BoundaryPointSet(atom *Walker)
|
|---|
| 23 | {
|
|---|
| 24 | node = Walker;
|
|---|
| 25 | LinesCount = 0;
|
|---|
| 26 | Nr = Walker->nr;
|
|---|
| 27 | }
|
|---|
| 28 | ;
|
|---|
| 29 |
|
|---|
| 30 | BoundaryPointSet::~BoundaryPointSet()
|
|---|
| 31 | {
|
|---|
| 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;
|
|---|
| 36 | }
|
|---|
| 37 | ;
|
|---|
| 38 |
|
|---|
| 39 | void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
|
|---|
| 40 | {
|
|---|
| 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 | else
|
|---|
| 48 | {
|
|---|
| 49 | lines.insert(LinePair(line->endpoints[0]->Nr, line));
|
|---|
| 50 | }
|
|---|
| 51 | LinesCount++;
|
|---|
| 52 | }
|
|---|
| 53 | ;
|
|---|
| 54 |
|
|---|
| 55 | ostream &
|
|---|
| 56 | operator <<(ostream &ost, BoundaryPointSet &a)
|
|---|
| 57 | {
|
|---|
| 58 | ost << "[" << a.Nr << "|" << a.node->Name << "]";
|
|---|
| 59 | return ost;
|
|---|
| 60 | }
|
|---|
| 61 | ;
|
|---|
| 62 |
|
|---|
| 63 | // ======================================== Lines on Boundary =================================
|
|---|
| 64 |
|
|---|
| 65 | BoundaryLineSet::BoundaryLineSet()
|
|---|
| 66 | {
|
|---|
| 67 | for (int i = 0; i < 2; i++)
|
|---|
| 68 | endpoints[i] = NULL;
|
|---|
| 69 | TrianglesCount = 0;
|
|---|
| 70 | Nr = -1;
|
|---|
| 71 | }
|
|---|
| 72 | ;
|
|---|
| 73 |
|
|---|
| 74 | BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
|
|---|
| 75 | {
|
|---|
| 76 | // set number
|
|---|
| 77 | Nr = number;
|
|---|
| 78 | // set endpoints in ascending order
|
|---|
| 79 | SetEndpointsOrdered(endpoints, Point[0], Point[1]);
|
|---|
| 80 | // add this line to the hash maps of both endpoints
|
|---|
| 81 | Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
|
|---|
| 82 | Point[1]->AddLine(this); //
|
|---|
| 83 | // clear triangles list
|
|---|
| 84 | TrianglesCount = 0;
|
|---|
| 85 | cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
|
|---|
| 86 | }
|
|---|
| 87 | ;
|
|---|
| 88 |
|
|---|
| 89 | BoundaryLineSet::~BoundaryLineSet()
|
|---|
| 90 | {
|
|---|
| 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 | } else
|
|---|
| 103 | cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
|
|---|
| 104 | } else
|
|---|
| 105 | 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;
|
|---|
| 109 | }
|
|---|
| 110 | ;
|
|---|
| 111 |
|
|---|
| 112 | void
|
|---|
| 113 | BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
|
|---|
| 114 | {
|
|---|
| 115 | cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
|
|---|
| 116 | << endl;
|
|---|
| 117 | triangles.insert(TrianglePair(triangle->Nr, triangle));
|
|---|
| 118 | TrianglesCount++;
|
|---|
| 119 | }
|
|---|
| 120 | ;
|
|---|
| 121 |
|
|---|
| 122 | ostream &
|
|---|
| 123 | operator <<(ostream &ost, BoundaryLineSet &a)
|
|---|
| 124 | {
|
|---|
| 125 | ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
|
|---|
| 126 | << a.endpoints[1]->node->Name << "]";
|
|---|
| 127 | return ost;
|
|---|
| 128 | }
|
|---|
| 129 | ;
|
|---|
| 130 |
|
|---|
| 131 | // ======================================== Triangles on Boundary =================================
|
|---|
| 132 |
|
|---|
| 133 |
|
|---|
| 134 | BoundaryTriangleSet::BoundaryTriangleSet()
|
|---|
| 135 | {
|
|---|
| 136 | for (int i = 0; i < 3; i++)
|
|---|
| 137 | {
|
|---|
| 138 | endpoints[i] = NULL;
|
|---|
| 139 | lines[i] = NULL;
|
|---|
| 140 | }
|
|---|
| 141 | Nr = -1;
|
|---|
| 142 | }
|
|---|
| 143 | ;
|
|---|
| 144 |
|
|---|
| 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;
|
|---|
| 184 | }
|
|---|
| 185 | ;
|
|---|
| 186 |
|
|---|
| 187 | BoundaryTriangleSet::~BoundaryTriangleSet()
|
|---|
| 188 | {
|
|---|
| 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 | } else
|
|---|
| 198 | cerr << "ERROR: This line " << i << " has already been free'd." << endl;
|
|---|
| 199 | } else
|
|---|
| 200 | cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
|
|---|
| 201 | }
|
|---|
| 202 | }
|
|---|
| 203 | ;
|
|---|
| 204 |
|
|---|
| 205 | void
|
|---|
| 206 | BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
|
|---|
| 207 | {
|
|---|
| 208 | // get normal vector
|
|---|
| 209 | 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.);
|
|---|
| 215 | }
|
|---|
| 216 | ;
|
|---|
| 217 |
|
|---|
| 218 | ostream &
|
|---|
| 219 | operator <<(ostream &ost, BoundaryTriangleSet &a)
|
|---|
| 220 | {
|
|---|
| 221 | ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
|
|---|
| 222 | << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
|
|---|
| 223 | return ost;
|
|---|
| 224 | }
|
|---|
| 225 | ;
|
|---|
| 226 |
|
|---|
| 227 | // ========================================== F U N C T I O N S =================================
|
|---|
| 228 |
|
|---|
| 229 | /** Finds the endpoint two lines are sharing.
|
|---|
| 230 | * \param *line1 first line
|
|---|
| 231 | * \param *line2 second line
|
|---|
| 232 | * \return point which is shared or NULL if none
|
|---|
| 233 | */
|
|---|
| 234 | class BoundaryPointSet *
|
|---|
| 235 | GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
|
|---|
| 236 | {
|
|---|
| 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 lines
|
|---|
| 244 | for (int j = 0; j < 2; j++)
|
|---|
| 245 | { // for both endpoints
|
|---|
| 246 | 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 endpoint
|
|---|
| 250 | node = OrderTest.first->second;
|
|---|
| 251 | cout << Verbose(5) << "Common endpoint of lines " << *line1
|
|---|
| 252 | << " and " << *line2 << " is: " << *node << "." << endl;
|
|---|
| 253 | j = 2;
|
|---|
| 254 | i = 2;
|
|---|
| 255 | break;
|
|---|
| 256 | }
|
|---|
| 257 | }
|
|---|
| 258 | return node;
|
|---|
| 259 | }
|
|---|
| 260 | ;
|
|---|
| 261 |
|
|---|
| 262 | /** Determines the boundary points of a cluster.
|
|---|
| 263 | * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
|
|---|
| 264 | * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
|
|---|
| 265 | * center and first and last point in the triple, it is thrown out.
|
|---|
| 266 | * \param *out output stream for debugging
|
|---|
| 267 | * \param *mol molecule structure representing the cluster
|
|---|
| 268 | */
|
|---|
| 269 | Boundaries *
|
|---|
| 270 | GetBoundaryPoints(ofstream *out, molecule *mol)
|
|---|
| 271 | {
|
|---|
| 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 axis
|
|---|
| 283 | 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 neighbours
|
|---|
| 300 | 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 side
|
|---|
| 308 | //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 | else
|
|---|
| 314 | angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
|
|---|
| 315 |
|
|---|
| 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 inserted
|
|---|
| 328 | }
|
|---|
| 329 | else
|
|---|
| 330 | { // same point exists, check first r, then distance of original vectors to center of gravity
|
|---|
| 331 | *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 slower
|
|---|
| 353 | BoundaryTestPair.first->second.second = Walker;
|
|---|
| 354 | *out << Verbose(2) << "Keeping new vector." << endl;
|
|---|
| 355 | }
|
|---|
| 356 | else
|
|---|
| 357 | {
|
|---|
| 358 | *out << Verbose(2) << "Keeping present vector." << endl;
|
|---|
| 359 | }
|
|---|
| 360 | }
|
|---|
| 361 | else
|
|---|
| 362 | {
|
|---|
| 363 | *out << Verbose(2) << "Keeping present vector." << endl;
|
|---|
| 364 | }
|
|---|
| 365 | }
|
|---|
| 366 | }
|
|---|
| 367 | // printing all inserted for debugging
|
|---|
| 368 | // {
|
|---|
| 369 | // *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
|
|---|
| 370 | // int i=0;
|
|---|
| 371 | // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
|
|---|
| 372 | // if (runner != BoundaryPoints[axis].begin())
|
|---|
| 373 | // *out << ", " << i << ": " << *runner->second.second;
|
|---|
| 374 | // else
|
|---|
| 375 | // *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 neighbours
|
|---|
| 381 | bool flag = false;
|
|---|
| 382 | do
|
|---|
| 383 | { // do as long as we still throw one out per round
|
|---|
| 384 | *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(); runner
|
|---|
| 391 | != BoundaryPoints[axis].end(); runner++)
|
|---|
| 392 | {
|
|---|
| 393 | // set neighbours correctly
|
|---|
| 394 | if (runner == BoundaryPoints[axis].begin())
|
|---|
| 395 | {
|
|---|
| 396 | left = BoundaryPoints[axis].end();
|
|---|
| 397 | }
|
|---|
| 398 | else
|
|---|
| 399 | {
|
|---|
| 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 distance
|
|---|
| 410 |
|
|---|
| 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 length
|
|---|
| 440 | double a = SideA.Norm();
|
|---|
| 441 | //double b = SideB.Norm();
|
|---|
| 442 | //double c = SideC.Norm();
|
|---|
| 443 | double h = SideH.Norm();
|
|---|
| 444 | // calculate the angles
|
|---|
| 445 | 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)) * (((alpha
|
|---|
| 450 | < 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 point
|
|---|
| 457 | //*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;
|
|---|
| 467 | }
|
|---|
| 468 | ;
|
|---|
| 469 |
|
|---|
| 470 | /** Determines greatest diameters of a cluster defined by its convex envelope.
|
|---|
| 471 | * Looks at lines parallel to one axis and where they intersect on the projected planes
|
|---|
| 472 | * \param *out output stream for debugging
|
|---|
| 473 | * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
|
|---|
| 474 | * \param *mol molecule structure representing the cluster
|
|---|
| 475 | * \param IsAngstroem whether we have angstroem or atomic units
|
|---|
| 476 | * \return NDIM array of the diameters
|
|---|
| 477 | */
|
|---|
| 478 | double *
|
|---|
| 479 | GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
|
|---|
| 480 | bool IsAngstroem)
|
|---|
| 481 | {
|
|---|
| 482 | // get points on boundary of NULL was given as parameter
|
|---|
| 483 | bool BoundaryFreeFlag = false;
|
|---|
| 484 | Boundaries *BoundaryPoints = BoundaryPtr;
|
|---|
| 485 | if (BoundaryPoints == NULL)
|
|---|
| 486 | {
|
|---|
| 487 | BoundaryFreeFlag = true;
|
|---|
| 488 | BoundaryPoints = GetBoundaryPoints(out, mol);
|
|---|
| 489 | }
|
|---|
| 490 | else
|
|---|
| 491 | {
|
|---|
| 492 | *out << Verbose(1) << "Using given boundary points set." << endl;
|
|---|
| 493 | }
|
|---|
| 494 | // determine biggest "diameter" of cluster for each axis
|
|---|
| 495 | 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 plane
|
|---|
| 504 | //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
|
|---|
| 505 | for (int j = 0; j < 2; j++)
|
|---|
| 506 | { // and for both axis on the current plane
|
|---|
| 507 | 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(); runner
|
|---|
| 511 | != 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 flips
|
|---|
| 515 | Neighbour = runner;
|
|---|
| 516 | Neighbour++;
|
|---|
| 517 | if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
|
|---|
| 518 | Neighbour = BoundaryPoints[axis].begin();
|
|---|
| 519 | DistanceVector.CopyVector(&runner->second.second->x);
|
|---|
| 520 | DistanceVector.SubtractVector(&Neighbour->second.second->x);
|
|---|
| 521 | do
|
|---|
| 522 | { // seek for neighbour pair where it flips
|
|---|
| 523 | OldComponent = DistanceVector.x[Othercomponent];
|
|---|
| 524 | Neighbour++;
|
|---|
| 525 | if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
|
|---|
| 526 | 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 flip
|
|---|
| 534 | if (runner != Neighbour)
|
|---|
| 535 | {
|
|---|
| 536 | OtherNeighbour = Neighbour;
|
|---|
| 537 | if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
|
|---|
| 538 | 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 OtherNeighbour
|
|---|
| 542 | 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 OtherNeighbour
|
|---|
| 547 | w1 = fabs(OtherVector.x[Othercomponent]);
|
|---|
| 548 | w2 = fabs(DistanceVector.x[Othercomponent]);
|
|---|
| 549 | tmp = fabs((w1 * DistanceVector.x[component] + w2
|
|---|
| 550 | * OtherVector.x[component]) / (w1 + w2));
|
|---|
| 551 | // mark if it has greater diameter
|
|---|
| 552 | //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
|
|---|
| 553 | GreatestDiameter[component] = (GreatestDiameter[component]
|
|---|
| 554 | > tmp) ? GreatestDiameter[component] : tmp;
|
|---|
| 555 | } //else
|
|---|
| 556 | //*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 lists
|
|---|
| 566 | if (BoundaryFreeFlag)
|
|---|
| 567 | delete[] (BoundaryPoints);
|
|---|
| 568 |
|
|---|
| 569 | return GreatestDiameter;
|
|---|
| 570 | }
|
|---|
| 571 | ;
|
|---|
| 572 |
|
|---|
| 573 | /** Creates the objects in a VRML file.
|
|---|
| 574 | * \param *out output stream for debugging
|
|---|
| 575 | * \param *vrmlfile output stream for tecplot data
|
|---|
| 576 | * \param *Tess Tesselation structure with constructed triangles
|
|---|
| 577 | * \param *mol molecule structure with atom positions
|
|---|
| 578 | */
|
|---|
| 579 | void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
|
|---|
| 580 | {
|
|---|
| 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 type
|
|---|
| 593 | 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 colour
|
|---|
| 596 | }
|
|---|
| 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 type
|
|---|
| 602 | 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 colour
|
|---|
| 608 | }
|
|---|
| 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 type
|
|---|
| 613 | for (i=0;i<3;i++) { // print each node
|
|---|
| 614 | for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates
|
|---|
| 615 | *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
|
|---|
| 616 | *vrmlfile << "\t";
|
|---|
| 617 | }
|
|---|
| 618 | *vrmlfile << "1. 0. 0." << endl; // red as colour
|
|---|
| 619 | *vrmlfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
|---|
| 620 | }
|
|---|
| 621 | } else {
|
|---|
| 622 | cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
|
|---|
| 623 | }
|
|---|
| 624 | delete(center);
|
|---|
| 625 | };
|
|---|
| 626 |
|
|---|
| 627 | /** Creates the objects in a raster3d file (renderable with a header.r3d).
|
|---|
| 628 | * \param *out output stream for debugging
|
|---|
| 629 | * \param *rasterfile output stream for tecplot data
|
|---|
| 630 | * \param *Tess Tesselation structure with constructed triangles
|
|---|
| 631 | * \param *mol molecule structure with atom positions
|
|---|
| 632 | */
|
|---|
| 633 | void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
|
|---|
| 634 | {
|
|---|
| 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 type
|
|---|
| 647 | 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 colour
|
|---|
| 650 | }
|
|---|
| 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 type
|
|---|
| 656 | 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 colour
|
|---|
| 662 | }
|
|---|
| 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.0 0 0\n";
|
|---|
| 666 | for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
|
|---|
| 667 | *rasterfile << "1" << endl << " "; // 1 is triangle type
|
|---|
| 668 | for (i=0;i<3;i++) { // print each node
|
|---|
| 669 | for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates
|
|---|
| 670 | *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
|
|---|
| 671 | *rasterfile << "\t";
|
|---|
| 672 | }
|
|---|
| 673 | *rasterfile << "1. 0. 0." << endl; // red as colour
|
|---|
| 674 | //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
|---|
| 675 | }
|
|---|
| 676 | *rasterfile << "9\n terminating special property\n";
|
|---|
| 677 | } else {
|
|---|
| 678 | cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
|
|---|
| 679 | }
|
|---|
| 680 | delete(center);
|
|---|
| 681 | };
|
|---|
| 682 |
|
|---|
| 683 | /** This function creates the tecplot file, displaying the tesselation of the hull.
|
|---|
| 684 | * \param *out output stream for debugging
|
|---|
| 685 | * \param *tecplot output stream for tecplot data
|
|---|
| 686 | * \param N arbitrary number to differentiate various zones in the tecplot format
|
|---|
| 687 | */
|
|---|
| 688 | void
|
|---|
| 689 | 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->TrianglesOnBoundaryCount
|
|---|
| 699 | << ", 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 coordinates
|
|---|
| 705 | *out << Verbose(2) << "The following triangles were created:";
|
|---|
| 706 | int Counter = 1;
|
|---|
| 707 | atom *Walker = NULL;
|
|---|
| 708 | for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target
|
|---|
| 709 | != 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 connectivity
|
|---|
| 718 | for (TriangleMap::iterator runner =
|
|---|
| 719 | TesselStruct->TrianglesOnBoundary.begin(); runner
|
|---|
| 720 | != 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 | }
|
|---|
| 732 | }
|
|---|
| 733 |
|
|---|
| 734 | /** Determines the volume of a cluster.
|
|---|
| 735 | * Determines first the convex envelope, then tesselates it and calculates its volume.
|
|---|
| 736 | * \param *out output stream for debugging
|
|---|
| 737 | * \param *filename filename prefix for output of vertex data
|
|---|
| 738 | * \param *configuration needed for path to store convex envelope file
|
|---|
| 739 | * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
|
|---|
| 740 | * \param *mol molecule structure representing the cluster
|
|---|
| 741 | * \return determined volume of the cluster in cubed config:GetIsAngstroem()
|
|---|
| 742 | */
|
|---|
| 743 | double
|
|---|
| 744 | 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 gravity
|
|---|
| 761 | *out << endl;
|
|---|
| 762 | Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
|
|---|
| 763 |
|
|---|
| 764 | // 2. translate all points into CoG
|
|---|
| 765 | *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 boundary
|
|---|
| 774 | if (BoundaryPoints == NULL)
|
|---|
| 775 | {
|
|---|
| 776 | BoundaryFreeFlag = true;
|
|---|
| 777 | BoundaryPoints = GetBoundaryPoints(out, mol);
|
|---|
| 778 | }
|
|---|
| 779 | else
|
|---|
| 780 | {
|
|---|
| 781 | *out << Verbose(1) << "Using given boundary points set." << endl;
|
|---|
| 782 | }
|
|---|
| 783 |
|
|---|
| 784 | // 4. fill the boundary point list
|
|---|
| 785 | for (int axis = 0; axis < NDIM; axis++)
|
|---|
| 786 | for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
|
|---|
| 787 | != BoundaryPoints[axis].end(); runner++)
|
|---|
| 788 | {
|
|---|
| 789 | TesselStruct->AddPoint(runner->second.second);
|
|---|
| 790 | }
|
|---|
| 791 |
|
|---|
| 792 | *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
|
|---|
| 793 | << " points on the convex boundary." << endl;
|
|---|
| 794 | // now we have the whole set of edge points in the BoundaryList
|
|---|
| 795 |
|
|---|
| 796 | // listing for debugging
|
|---|
| 797 | // *out << Verbose(1) << "Listing PointsOnBoundary:";
|
|---|
| 798 | // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
|
|---|
| 799 | // *out << " " << *runner->second;
|
|---|
| 800 | // }
|
|---|
| 801 | // *out << endl;
|
|---|
| 802 |
|
|---|
| 803 | // 5a. guess starting triangle
|
|---|
| 804 | 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->TrianglesOnBoundaryCount
|
|---|
| 810 | << " triangles with " << TesselStruct->LinesOnBoundaryCount
|
|---|
| 811 | << " 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 volumes
|
|---|
| 815 | *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(); runner
|
|---|
| 819 | != TesselStruct->TrianglesOnBoundary.end(); runner++)
|
|---|
| 820 | { // go through every triangle, calculate volume of its pyramid with CoG as peak
|
|---|
| 821 | 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 triangle
|
|---|
| 832 | 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 triangle
|
|---|
| 837 | 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 CoG
|
|---|
| 849 | *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 file
|
|---|
| 860 | 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 lists
|
|---|
| 868 | if (BoundaryFreeFlag)
|
|---|
| 869 | delete[] (BoundaryPoints);
|
|---|
| 870 |
|
|---|
| 871 | return volume;
|
|---|
| 872 | }
|
|---|
| 873 | ;
|
|---|
| 874 |
|
|---|
| 875 | /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
|
|---|
| 876 | * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
|
|---|
| 877 | * \param *out output stream for debugging
|
|---|
| 878 | * \param *configuration needed for path to store convex envelope file
|
|---|
| 879 | * \param *mol molecule structure representing the cluster
|
|---|
| 880 | * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
|
|---|
| 881 | * \param celldensity desired average density in final cell
|
|---|
| 882 | */
|
|---|
| 883 | void
|
|---|
| 884 | PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
|
|---|
| 885 | double ClusterVolume, double celldensity)
|
|---|
| 886 | {
|
|---|
| 887 | // transform to PAS
|
|---|
| 888 | mol->PrincipalAxisSystem(out, true);
|
|---|
| 889 |
|
|---|
| 890 | // some preparations beforehand
|
|---|
| 891 | 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 | else
|
|---|
| 898 | 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 masses
|
|---|
| 909 | 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 polynomial
|
|---|
| 924 | *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
|
|---|
| 925 | << endl;
|
|---|
| 926 | double cellvolume;
|
|---|
| 927 | if (IsAngstroem)
|
|---|
| 928 | cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
|
|---|
| 929 | / clustervolume)) / (celldensity - 1);
|
|---|
| 930 | else
|
|---|
| 931 | cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
|
|---|
| 932 | / clustervolume)) / (celldensity - 1);
|
|---|
| 933 | *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
|
|---|
| 934 | << " 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 | else
|
|---|
| 956 | {
|
|---|
| 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 return
|
|---|
| 967 | *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
|
|---|
| 968 | << " ." << endl;
|
|---|
| 969 | else
|
|---|
| 970 | {
|
|---|
| 971 | *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
|
|---|
| 972 | << " and " << x1 << " and " << x2 << " ." << endl;
|
|---|
| 973 | x0 = x2; // sorted in ascending order
|
|---|
| 974 | }
|
|---|
| 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 dimensions
|
|---|
| 984 | *out << Verbose(0) << "Translating to box with these boundaries." << endl;
|
|---|
| 985 | mol->CenterInBox((ofstream *) &cout, &BoxLengths);
|
|---|
| 986 | }
|
|---|
| 987 | // update Box of atoms by boundary
|
|---|
| 988 | 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;
|
|---|
| 993 | }
|
|---|
| 994 | ;
|
|---|
| 995 |
|
|---|
| 996 | // =========================================================== class TESSELATION ===========================================
|
|---|
| 997 |
|
|---|
| 998 | /** Constructor of class Tesselation.
|
|---|
| 999 | */
|
|---|
| 1000 | Tesselation::Tesselation()
|
|---|
| 1001 | {
|
|---|
| 1002 | PointsOnBoundaryCount = 0;
|
|---|
| 1003 | LinesOnBoundaryCount = 0;
|
|---|
| 1004 | TrianglesOnBoundaryCount = 0;
|
|---|
| 1005 | TriangleFilesWritten = 0;
|
|---|
| 1006 | }
|
|---|
| 1007 | ;
|
|---|
| 1008 |
|
|---|
| 1009 | /** Constructor of class Tesselation.
|
|---|
| 1010 | * We have to free all points, lines and triangles.
|
|---|
| 1011 | */
|
|---|
| 1012 | Tesselation::~Tesselation()
|
|---|
| 1013 | {
|
|---|
| 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 | } else
|
|---|
| 1020 | cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
|
|---|
| 1021 | }
|
|---|
| 1022 | }
|
|---|
| 1023 | ;
|
|---|
| 1024 |
|
|---|
| 1025 | /** Gueses first starting triangle of the convex envelope.
|
|---|
| 1026 | * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
|
|---|
| 1027 | * \param *out output stream for debugging
|
|---|
| 1028 | * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
|
|---|
| 1029 | */
|
|---|
| 1030 | void
|
|---|
| 1031 | Tesselation::GuessStartingTriangle(ofstream *out)
|
|---|
| 1032 | {
|
|---|
| 1033 | // 4b. create a starting triangle
|
|---|
| 1034 | // 4b1. create all distances
|
|---|
| 1035 | DistanceMultiMap DistanceMMap;
|
|---|
| 1036 | double distance, tmp;
|
|---|
| 1037 | Vector PlaneVector, TrialVector;
|
|---|
| 1038 | PointMap::iterator A, B, C; // three nodes of the first triangle
|
|---|
| 1039 | A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
|
|---|
| 1040 |
|
|---|
| 1041 | // with A chosen, take each pair B,C and sort
|
|---|
| 1042 | 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 distances
|
|---|
| 1064 | // *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 triangle
|
|---|
| 1070 | // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
|
|---|
| 1071 | 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 candidate
|
|---|
| 1075 | // 2. next, we have to check whether all points reside on only one side of the triangle
|
|---|
| 1076 | // 3. construct plane vector
|
|---|
| 1077 | 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 points
|
|---|
| 1084 | 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) || (checker
|
|---|
| 1090 | == baseline->second.second))
|
|---|
| 1091 | continue;
|
|---|
| 1092 | // 4a. project onto plane vector
|
|---|
| 1093 | 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 ok
|
|---|
| 1097 | continue;
|
|---|
| 1098 | *out << Verbose(3) << "Projection of " << checker->second->node->Name
|
|---|
| 1099 | << " 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. loop
|
|---|
| 1105 | *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 | else
|
|---|
| 1114 | { // note the sign for later
|
|---|
| 1115 | *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 node
|
|---|
| 1124 | 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)) && (tmp
|
|---|
| 1128 | < 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)) && (tmp
|
|---|
| 1135 | < 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)) && (tmp
|
|---|
| 1142 | < 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. loop
|
|---|
| 1146 | if (innerpoint == 3)
|
|---|
| 1147 | break;
|
|---|
| 1148 | }
|
|---|
| 1149 | // 5. come this far, all on same side? Then break 1. loop and construct triangle
|
|---|
| 1150 | 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 triangle
|
|---|
| 1169 | 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 | else
|
|---|
| 1181 | {
|
|---|
| 1182 | *out << Verbose(1) << "No starting triangle found." << endl;
|
|---|
| 1183 | exit(255);
|
|---|
| 1184 | }
|
|---|
| 1185 | }
|
|---|
| 1186 | ;
|
|---|
| 1187 |
|
|---|
| 1188 | /** Tesselates the convex envelope of a cluster from a single starting triangle.
|
|---|
| 1189 | * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
|
|---|
| 1190 | * 2 triangles. Hence, we go through all current lines:
|
|---|
| 1191 | * -# if the lines contains to only one triangle
|
|---|
| 1192 | * -# 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 vectors
|
|---|
| 1194 | * -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
|
|---|
| 1195 | * 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)
|
|---|
| 1197 | * \param *out output stream for debugging
|
|---|
| 1198 | * \param *configuration for IsAngstroem
|
|---|
| 1199 | * \param *mol the cluster as a molecule structure
|
|---|
| 1200 | */
|
|---|
| 1201 | void
|
|---|
| 1202 | 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 | do
|
|---|
| 1213 | {
|
|---|
| 1214 | flag = false;
|
|---|
| 1215 | for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
|
|---|
| 1216 | != 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 far
|
|---|
| 1224 | // get peak point with respect to this base line's only triangle
|
|---|
| 1225 | 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 triangle
|
|---|
| 1231 | BTS->GetNormalVector(NormalVector);
|
|---|
| 1232 | *out << Verbose(4) << "NormalVector of base triangle is ";
|
|---|
| 1233 | NormalVector.Output(out);
|
|---|
| 1234 | *out << endl;
|
|---|
| 1235 | // offset to center of triangle
|
|---|
| 1236 | 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 baseline
|
|---|
| 1252 | 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(); target
|
|---|
| 1258 | != PointsOnBoundary.end(); target++)
|
|---|
| 1259 | if ((target->second != baseline->second->endpoints[0])
|
|---|
| 1260 | && (target->second != baseline->second->endpoints[1]))
|
|---|
| 1261 | { // don't take the same endpoints
|
|---|
| 1262 | *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 line
|
|---|
| 1271 | VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
|
|---|
| 1272 | 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 " << TempAngle
|
|---|
| 1279 | << ", bad direction!" << endl;
|
|---|
| 1280 | continue;
|
|---|
| 1281 | }
|
|---|
| 1282 | else
|
|---|
| 1283 | *out << Verbose(4)
|
|---|
| 1284 | << "Angle between propagation direction and base line to "
|
|---|
| 1285 | << *(target->second) << " is " << TempAngle
|
|---|
| 1286 | << ", 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 | // else
|
|---|
| 1294 | // *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 | // else
|
|---|
| 1298 | // *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->TrianglesCount
|
|---|
| 1310 | << " 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->TrianglesCount
|
|---|
| 1324 | << " 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 found
|
|---|
| 1340 | if (continueflag)
|
|---|
| 1341 | { // create virtually this triangle, get its normal vector, calculate angle
|
|---|
| 1342 | 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 inward
|
|---|
| 1348 | if (baseline->second->endpoints[0]->node->x.Projection(
|
|---|
| 1349 | &VirtualNormalVector) > 0)
|
|---|
| 1350 | VirtualNormalVector.Scale(-1.);
|
|---|
| 1351 | // calculate angle
|
|---|
| 1352 | 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 winner
|
|---|
| 1358 | 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 triangle
|
|---|
| 1364 | if (winner != PointsOnBoundary.end())
|
|---|
| 1365 | {
|
|---|
| 1366 | *out << Verbose(2) << "Winning target point is "
|
|---|
| 1367 | << *(winner->second) << " with angle " << SmallestAngle
|
|---|
| 1368 | << "." << endl;
|
|---|
| 1369 | // create the lins of not yet present
|
|---|
| 1370 | 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 | { // create
|
|---|
| 1379 | 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 | else
|
|---|
| 1388 | BLS[1] = LineChecker[0]->second;
|
|---|
| 1389 | if (LineChecker[1]
|
|---|
| 1390 | == baseline->second->endpoints[1]->lines.end())
|
|---|
| 1391 | { // create
|
|---|
| 1392 | 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 | else
|
|---|
| 1401 | BLS[2] = LineChecker[1]->second;
|
|---|
| 1402 | BTS = new class BoundaryTriangleSet(BLS,
|
|---|
| 1403 | TrianglesOnBoundaryCount);
|
|---|
| 1404 | TrianglesOnBoundary.insert(TrianglePair(
|
|---|
| 1405 | TrianglesOnBoundaryCount, BTS));
|
|---|
| 1406 | TrianglesOnBoundaryCount++;
|
|---|
| 1407 | }
|
|---|
| 1408 | else
|
|---|
| 1409 | {
|
|---|
| 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 continue
|
|---|
| 1416 | }
|
|---|
| 1417 | else
|
|---|
| 1418 | *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
|
|---|
| 1419 | << " has a triangle count of "
|
|---|
| 1420 | << baseline->second->TrianglesCount << "." << endl;
|
|---|
| 1421 | }
|
|---|
| 1422 | while (flag);
|
|---|
| 1423 |
|
|---|
| 1424 | }
|
|---|
| 1425 | ;
|
|---|
| 1426 |
|
|---|
| 1427 | /** Adds an atom to the tesselation::PointsOnBoundary list.
|
|---|
| 1428 | * \param *Walker atom to add
|
|---|
| 1429 | */
|
|---|
| 1430 | void
|
|---|
| 1431 | Tesselation::AddPoint(atom *Walker)
|
|---|
| 1432 | {
|
|---|
| 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 counter
|
|---|
| 1437 | PointsOnBoundaryCount++;
|
|---|
| 1438 | }
|
|---|
| 1439 | ;
|
|---|
| 1440 |
|
|---|
| 1441 | /** Adds point to Tesselation::PointsOnBoundary if not yet present.
|
|---|
| 1442 | * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
|
|---|
| 1443 | * @param Candidate point to add
|
|---|
| 1444 | * @param n index for this point in Tesselation::TPS array
|
|---|
| 1445 | */
|
|---|
| 1446 | void
|
|---|
| 1447 | Tesselation::AddTrianglePoint(atom* Candidate, int n)
|
|---|
| 1448 | {
|
|---|
| 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 counter
|
|---|
| 1453 | {
|
|---|
| 1454 | PointsOnBoundaryCount++;
|
|---|
| 1455 | }
|
|---|
| 1456 | else
|
|---|
| 1457 | {
|
|---|
| 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 | }
|
|---|
| 1463 | }
|
|---|
| 1464 | ;
|
|---|
| 1465 |
|
|---|
| 1466 | /** Function tries to add line from current Points in BPS to BoundaryLineSet.
|
|---|
| 1467 | * If succesful it raises the line count and inserts the new line into the BLS,
|
|---|
| 1468 | * if unsuccesful, it writes the line which had been present into the BLS, deleting the new constructed one.
|
|---|
| 1469 | * @param *a first endpoint
|
|---|
| 1470 | * @param *b second endpoint
|
|---|
| 1471 | * @param n index of Tesselation::BLS giving the line with both endpoints
|
|---|
| 1472 | */
|
|---|
| 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 | }
|
|---|
| 1512 | }
|
|---|
| 1513 | ;
|
|---|
| 1514 |
|
|---|
| 1515 | /** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
|
|---|
| 1516 | * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
|
|---|
| 1517 | */
|
|---|
| 1518 | void
|
|---|
| 1519 | Tesselation::AddTriangleToLines()
|
|---|
| 1520 | {
|
|---|
| 1521 |
|
|---|
| 1522 | cout << Verbose(1) << "Adding triangle to its lines" << endl;
|
|---|
| 1523 | TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
|
|---|
| 1524 | TrianglesOnBoundaryCount++;
|
|---|
| 1525 |
|
|---|
| 1526 | /*
|
|---|
| 1527 | * this is apparently done when constructing triangle
|
|---|
| 1528 |
|
|---|
| 1529 | for (i=0; i<3; i++)
|
|---|
| 1530 | {
|
|---|
| 1531 | BLS[i]->AddTriangle(BTS);
|
|---|
| 1532 | }
|
|---|
| 1533 | */
|
|---|
| 1534 | }
|
|---|
| 1535 | ;
|
|---|
| 1536 |
|
|---|
| 1537 |
|
|---|
| 1538 | double det_get(gsl_matrix *A, int inPlace) {
|
|---|
| 1539 | /*
|
|---|
| 1540 | inPlace = 1 => A is replaced with the LU decomposed copy.
|
|---|
| 1541 | inPlace = 0 => A is retained, and a copy is used for LU.
|
|---|
| 1542 | */
|
|---|
| 1543 |
|
|---|
| 1544 | double det;
|
|---|
| 1545 | int signum;
|
|---|
| 1546 | gsl_permutation *p = gsl_permutation_alloc(A->size1);
|
|---|
| 1547 | gsl_matrix *tmpA;
|
|---|
| 1548 |
|
|---|
| 1549 | if (inPlace)
|
|---|
| 1550 | tmpA = A;
|
|---|
| 1551 | else {
|
|---|
| 1552 | gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
|
|---|
| 1553 | gsl_matrix_memcpy(tmpA , A);
|
|---|
| 1554 | }
|
|---|
| 1555 |
|
|---|
| 1556 |
|
|---|
| 1557 | gsl_linalg_LU_decomp(tmpA , p , &signum);
|
|---|
| 1558 | det = gsl_linalg_LU_det(tmpA , signum);
|
|---|
| 1559 | gsl_permutation_free(p);
|
|---|
| 1560 | if (! inPlace)
|
|---|
| 1561 | gsl_matrix_free(tmpA);
|
|---|
| 1562 |
|
|---|
| 1563 | return det;
|
|---|
| 1564 | };
|
|---|
| 1565 |
|
|---|
| 1566 | void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
|
|---|
| 1567 | {
|
|---|
| 1568 | gsl_matrix *A = gsl_matrix_calloc(3,3);
|
|---|
| 1569 | double m11, m12, m13, m14;
|
|---|
| 1570 |
|
|---|
| 1571 | for(int i=0;i<3;i++) {
|
|---|
| 1572 | gsl_matrix_set(A, i, 0, a.x[i]);
|
|---|
| 1573 | gsl_matrix_set(A, i, 1, b.x[i]);
|
|---|
| 1574 | gsl_matrix_set(A, i, 2, c.x[i]);
|
|---|
| 1575 | }
|
|---|
| 1576 | m11 = det_get(A, 1);
|
|---|
| 1577 |
|
|---|
| 1578 | for(int i=0;i<3;i++) {
|
|---|
| 1579 | gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
|
|---|
| 1580 | gsl_matrix_set(A, i, 1, b.x[i]);
|
|---|
| 1581 | gsl_matrix_set(A, i, 2, c.x[i]);
|
|---|
| 1582 | }
|
|---|
| 1583 | m12 = det_get(A, 1);
|
|---|
| 1584 |
|
|---|
| 1585 | for(int i=0;i<3;i++) {
|
|---|
| 1586 | gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
|
|---|
| 1587 | gsl_matrix_set(A, i, 1, a.x[i]);
|
|---|
| 1588 | gsl_matrix_set(A, i, 2, c.x[i]);
|
|---|
| 1589 | }
|
|---|
| 1590 | m13 = det_get(A, 1);
|
|---|
| 1591 |
|
|---|
| 1592 | for(int i=0;i<3;i++) {
|
|---|
| 1593 | gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
|
|---|
| 1594 | gsl_matrix_set(A, i, 1, a.x[i]);
|
|---|
| 1595 | gsl_matrix_set(A, i, 2, b.x[i]);
|
|---|
| 1596 | }
|
|---|
| 1597 | m14 = det_get(A, 1);
|
|---|
| 1598 |
|
|---|
| 1599 | if (fabs(m11) < MYEPSILON)
|
|---|
| 1600 | cerr << "ERROR: three points are colinear." << endl;
|
|---|
| 1601 |
|
|---|
| 1602 | center->x[0] = 0.5 * m12/ m11;
|
|---|
| 1603 | center->x[1] = -0.5 * m13/ m11;
|
|---|
| 1604 | center->x[2] = 0.5 * m14/ m11;
|
|---|
| 1605 |
|
|---|
| 1606 | if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
|
|---|
| 1607 | cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
|
|---|
| 1608 |
|
|---|
| 1609 | gsl_matrix_free(A);
|
|---|
| 1610 | };
|
|---|
| 1611 |
|
|---|
| 1612 |
|
|---|
| 1613 |
|
|---|
| 1614 | /**
|
|---|
| 1615 | * Function returns center of sphere with RADIUS, which rests on points a, b, c
|
|---|
| 1616 | * @param Center this vector will be used for return
|
|---|
| 1617 | * @param a vector first point of triangle
|
|---|
| 1618 | * @param b vector second point of triangle
|
|---|
| 1619 | * @param c vector third point of triangle
|
|---|
| 1620 | * @param Direction vector indicates up/down
|
|---|
| 1621 | * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
|
|---|
| 1622 | * @param Halfplaneindicator double indicates whether Direction is up or down
|
|---|
| 1623 | * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
|
|---|
| 1624 | * @param alpha double angle at a
|
|---|
| 1625 | * @param beta double, angle at b
|
|---|
| 1626 | * @param gamma, double, angle at c
|
|---|
| 1627 | * @param Radius, double
|
|---|
| 1628 | * @param Umkreisradius double radius of circumscribing circle
|
|---|
| 1629 | */
|
|---|
| 1630 | void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
|
|---|
| 1631 | double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
|
|---|
| 1632 | {
|
|---|
| 1633 | Vector TempNormal, helper;
|
|---|
| 1634 | double Restradius;
|
|---|
| 1635 | Vector OtherCenter;
|
|---|
| 1636 | cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
|
|---|
| 1637 | Center->Zero();
|
|---|
| 1638 | helper.CopyVector(&a);
|
|---|
| 1639 | helper.Scale(sin(2.*alpha));
|
|---|
| 1640 | Center->AddVector(&helper);
|
|---|
| 1641 | helper.CopyVector(&b);
|
|---|
| 1642 | helper.Scale(sin(2.*beta));
|
|---|
| 1643 | Center->AddVector(&helper);
|
|---|
| 1644 | helper.CopyVector(&c);
|
|---|
| 1645 | helper.Scale(sin(2.*gamma));
|
|---|
| 1646 | Center->AddVector(&helper);
|
|---|
| 1647 | //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
|
|---|
| 1648 | Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
|
|---|
| 1649 | NewUmkreismittelpunkt->CopyVector(Center);
|
|---|
| 1650 | cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
|
|---|
| 1651 | // Here we calculated center of circumscribing circle, using barycentric coordinates
|
|---|
| 1652 | cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
|
|---|
| 1653 |
|
|---|
| 1654 | TempNormal.CopyVector(&a);
|
|---|
| 1655 | TempNormal.SubtractVector(&b);
|
|---|
| 1656 | helper.CopyVector(&a);
|
|---|
| 1657 | helper.SubtractVector(&c);
|
|---|
| 1658 | TempNormal.VectorProduct(&helper);
|
|---|
| 1659 | if (fabs(HalfplaneIndicator) < MYEPSILON)
|
|---|
| 1660 | {
|
|---|
| 1661 | if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
|
|---|
| 1662 | {
|
|---|
| 1663 | TempNormal.Scale(-1);
|
|---|
| 1664 | }
|
|---|
| 1665 | }
|
|---|
| 1666 | else
|
|---|
| 1667 | {
|
|---|
| 1668 | if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
|
|---|
| 1669 | {
|
|---|
| 1670 | TempNormal.Scale(-1);
|
|---|
| 1671 | }
|
|---|
| 1672 | }
|
|---|
| 1673 |
|
|---|
| 1674 | TempNormal.Normalize();
|
|---|
| 1675 | Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
|
|---|
| 1676 | cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
|
|---|
| 1677 | TempNormal.Scale(Restradius);
|
|---|
| 1678 | cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
|
|---|
| 1679 |
|
|---|
| 1680 | Center->AddVector(&TempNormal);
|
|---|
| 1681 | cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
|
|---|
| 1682 | get_sphere(&OtherCenter, a, b, c, RADIUS);
|
|---|
| 1683 | cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
|
|---|
| 1684 | cout << Verbose(3) << "End of Get_center_of_sphere.\n";
|
|---|
| 1685 | };
|
|---|
| 1686 |
|
|---|
| 1687 | /** This recursive function finds a third point, to form a triangle with two given ones.
|
|---|
| 1688 | * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
|
|---|
| 1689 | * supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
|
|---|
| 1690 | * upon which we operate.
|
|---|
| 1691 | * If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
|
|---|
| 1692 | * direction and angle into Storage.
|
|---|
| 1693 | * We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
|
|---|
| 1694 | * with all neighbours of the candidate.
|
|---|
| 1695 | * @param a first point
|
|---|
| 1696 | * @param b second point
|
|---|
| 1697 | * *param c atom old third point of old triangle
|
|---|
| 1698 | * @param Candidate base point along whose bonds to start looking from
|
|---|
| 1699 | * @param Parent point to avoid during search as its wrong direction
|
|---|
| 1700 | * @param RecursionLevel contains current recursion depth
|
|---|
| 1701 | * @param Chord baseline vector of first and second point
|
|---|
| 1702 | * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
|
|---|
| 1703 | * @param OldNormal normal of the triangle which the baseline belongs to
|
|---|
| 1704 | * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere
|
|---|
| 1705 | * @param Opt_Candidate candidate reference to return
|
|---|
| 1706 | * @param Storage array containing two angles of current Opt_Candidate
|
|---|
| 1707 | * @param RADIUS radius of ball
|
|---|
| 1708 | * @param mol molecule structure with atoms and bonds
|
|---|
| 1709 | */
|
|---|
| 1710 | void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
|
|---|
| 1711 | int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
|
|---|
| 1712 | atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
|
|---|
| 1713 | {
|
|---|
| 1714 | cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
|
|---|
| 1715 | cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
|
|---|
| 1716 | cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
|
|---|
| 1717 | cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
|
|---|
| 1718 | cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
|
|---|
| 1719 | cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
|
|---|
| 1720 | /* OldNormal is normal vector on the old triangle
|
|---|
| 1721 | * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
|
|---|
| 1722 | * Chord points from b to a!!!
|
|---|
| 1723 | */
|
|---|
| 1724 | Vector dif_a; //Vector from a to candidate
|
|---|
| 1725 | Vector dif_b; //Vector from b to candidate
|
|---|
| 1726 | Vector AngleCheck;
|
|---|
| 1727 | Vector TempNormal, Umkreismittelpunkt;
|
|---|
| 1728 | Vector Mittelpunkt;
|
|---|
| 1729 |
|
|---|
| 1730 | double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
|
|---|
| 1731 | double BallAngle, AlternativeSign;
|
|---|
| 1732 | atom *Walker; // variable atom point
|
|---|
| 1733 |
|
|---|
| 1734 | Vector NewUmkreismittelpunkt;
|
|---|
| 1735 |
|
|---|
| 1736 | if (a != Candidate and b != Candidate and c != Candidate) {
|
|---|
| 1737 | cout << Verbose(3) << "We have a unique candidate!" << endl;
|
|---|
| 1738 | dif_a.CopyVector(&(a->x));
|
|---|
| 1739 | dif_a.SubtractVector(&(Candidate->x));
|
|---|
| 1740 | dif_b.CopyVector(&(b->x));
|
|---|
| 1741 | dif_b.SubtractVector(&(Candidate->x));
|
|---|
| 1742 | AngleCheck.CopyVector(&(Candidate->x));
|
|---|
| 1743 | AngleCheck.SubtractVector(&(a->x));
|
|---|
| 1744 | AngleCheck.ProjectOntoPlane(Chord);
|
|---|
| 1745 |
|
|---|
| 1746 | SideA = dif_b.Norm();
|
|---|
| 1747 | SideB = dif_a.Norm();
|
|---|
| 1748 | SideC = Chord->Norm();
|
|---|
| 1749 | //Chord->Scale(-1);
|
|---|
| 1750 |
|
|---|
| 1751 | alpha = Chord->Angle(&dif_a);
|
|---|
| 1752 | beta = M_PI - Chord->Angle(&dif_b);
|
|---|
| 1753 | gamma = dif_a.Angle(&dif_b);
|
|---|
| 1754 |
|
|---|
| 1755 | 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;
|
|---|
| 1756 |
|
|---|
| 1757 | if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
|
|---|
| 1758 | cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
|
|---|
| 1759 | 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;
|
|---|
| 1760 | }
|
|---|
| 1761 |
|
|---|
| 1762 | if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
|
|---|
| 1763 | Umkreisradius = SideA / 2.0 / sin(alpha);
|
|---|
| 1764 | //cout << Umkreisradius << endl;
|
|---|
| 1765 | //cout << SideB / 2.0 / sin(beta) << endl;
|
|---|
| 1766 | //cout << SideC / 2.0 / sin(gamma) << endl;
|
|---|
| 1767 |
|
|---|
| 1768 | if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
|
|---|
| 1769 | cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
|
|---|
| 1770 | cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
|
|---|
| 1771 | sign = AngleCheck.ScalarProduct(direction1);
|
|---|
| 1772 | if (fabs(sign)<MYEPSILON) {
|
|---|
| 1773 | if (AngleCheck.ScalarProduct(OldNormal)<0) {
|
|---|
| 1774 | sign =0;
|
|---|
| 1775 | AlternativeSign=1;
|
|---|
| 1776 | } else {
|
|---|
| 1777 | sign =0;
|
|---|
| 1778 | AlternativeSign=-1;
|
|---|
| 1779 | }
|
|---|
| 1780 | } else {
|
|---|
| 1781 | sign /= fabs(sign);
|
|---|
| 1782 | }
|
|---|
| 1783 | if (sign >= 0) {
|
|---|
| 1784 | cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
|
|---|
| 1785 | Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
|
|---|
| 1786 | Mittelpunkt.SubtractVector(&ReferencePoint);
|
|---|
| 1787 | cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
|
|---|
| 1788 | BallAngle = Mittelpunkt.Angle(OldNormal);
|
|---|
| 1789 | cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
|
|---|
| 1790 |
|
|---|
| 1791 | //cout << "direction1 is " << *direction1 << "." << endl;
|
|---|
| 1792 | //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
|
|---|
| 1793 | //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
|
|---|
| 1794 |
|
|---|
| 1795 | NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
|
|---|
| 1796 |
|
|---|
| 1797 | if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
|
|---|
| 1798 | if (Storage[0]< -1.5) { // first Candidate at all
|
|---|
| 1799 | if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
|
|---|
| 1800 | cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
|
|---|
| 1801 | Opt_Candidate = Candidate;
|
|---|
| 1802 | Storage[0] = sign;
|
|---|
| 1803 | Storage[1] = AlternativeSign;
|
|---|
| 1804 | Storage[2] = BallAngle;
|
|---|
| 1805 | cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
|
|---|
| 1806 | } else
|
|---|
| 1807 | cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
|
|---|
| 1808 | } else {
|
|---|
| 1809 | if ( Storage[2] > BallAngle) {
|
|---|
| 1810 | if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
|
|---|
| 1811 | cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
|
|---|
| 1812 | Opt_Candidate = Candidate;
|
|---|
| 1813 | Storage[0] = sign;
|
|---|
| 1814 | Storage[1] = AlternativeSign;
|
|---|
| 1815 | Storage[2] = BallAngle;
|
|---|
| 1816 | cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
|
|---|
| 1817 | } else
|
|---|
| 1818 | cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
|
|---|
| 1819 | } else {
|
|---|
| 1820 | if (DEBUG) {
|
|---|
| 1821 | cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
|
|---|
| 1822 | }
|
|---|
| 1823 | }
|
|---|
| 1824 | }
|
|---|
| 1825 | } else {
|
|---|
| 1826 | if (DEBUG) {
|
|---|
| 1827 | cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
|
|---|
| 1828 | }
|
|---|
| 1829 | }
|
|---|
| 1830 | } else {
|
|---|
| 1831 | if (DEBUG) {
|
|---|
| 1832 | cout << Verbose(3) << *Candidate << " is not in search direction." << endl;
|
|---|
| 1833 | }
|
|---|
| 1834 | }
|
|---|
| 1835 | } else {
|
|---|
| 1836 | if (DEBUG) {
|
|---|
| 1837 | cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
|
|---|
| 1838 | }
|
|---|
| 1839 | }
|
|---|
| 1840 | } else {
|
|---|
| 1841 | if (DEBUG) {
|
|---|
| 1842 | cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
|
|---|
| 1843 | }
|
|---|
| 1844 | }
|
|---|
| 1845 | } else {
|
|---|
| 1846 | if (DEBUG) {
|
|---|
| 1847 | cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
|
|---|
| 1848 | }
|
|---|
| 1849 | }
|
|---|
| 1850 |
|
|---|
| 1851 | if (RecursionLevel < 5) { // Seven is the recursion level threshold.
|
|---|
| 1852 | for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
|
|---|
| 1853 | Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
|
|---|
| 1854 | if (Walker == Parent) { // don't go back the same bond
|
|---|
| 1855 | continue;
|
|---|
| 1856 | } else {
|
|---|
| 1857 | 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
|
|---|
| 1858 | }
|
|---|
| 1859 | }
|
|---|
| 1860 | }
|
|---|
| 1861 | cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
|
|---|
| 1862 | }
|
|---|
| 1863 | ;
|
|---|
| 1864 |
|
|---|
| 1865 |
|
|---|
| 1866 | /** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
|
|---|
| 1867 | * \param *Center new center on return
|
|---|
| 1868 | * \param *a first point
|
|---|
| 1869 | * \param *b second point
|
|---|
| 1870 | * \param *c third point
|
|---|
| 1871 | */
|
|---|
| 1872 | void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
|
|---|
| 1873 | {
|
|---|
| 1874 | Vector helper;
|
|---|
| 1875 | double alpha, beta, gamma;
|
|---|
| 1876 | Vector SideA, SideB, SideC;
|
|---|
| 1877 | SideA.CopyVector(b);
|
|---|
| 1878 | SideA.SubtractVector(c);
|
|---|
| 1879 | SideB.CopyVector(c);
|
|---|
| 1880 | SideB.SubtractVector(a);
|
|---|
| 1881 | SideC.CopyVector(a);
|
|---|
| 1882 | SideC.SubtractVector(b);
|
|---|
| 1883 | alpha = M_PI - SideB.Angle(&SideC);
|
|---|
| 1884 | beta = M_PI - SideC.Angle(&SideA);
|
|---|
| 1885 | gamma = M_PI - SideA.Angle(&SideB);
|
|---|
| 1886 | cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
|
|---|
| 1887 | if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
|
|---|
| 1888 | cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
|
|---|
| 1889 |
|
|---|
| 1890 | Center->Zero();
|
|---|
| 1891 | helper.CopyVector(a);
|
|---|
| 1892 | helper.Scale(sin(2.*alpha));
|
|---|
| 1893 | Center->AddVector(&helper);
|
|---|
| 1894 | helper.CopyVector(b);
|
|---|
| 1895 | helper.Scale(sin(2.*beta));
|
|---|
| 1896 | Center->AddVector(&helper);
|
|---|
| 1897 | helper.CopyVector(c);
|
|---|
| 1898 | helper.Scale(sin(2.*gamma));
|
|---|
| 1899 | Center->AddVector(&helper);
|
|---|
| 1900 | Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
|
|---|
| 1901 | };
|
|---|
| 1902 |
|
|---|
| 1903 | /** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
|
|---|
| 1904 | * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
|
|---|
| 1905 | * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
|
|---|
| 1906 | * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
|
|---|
| 1907 | * \param CircleCenter Center of the parameter circle
|
|---|
| 1908 | * \param CirclePlaneNormal normal vector to plane of the parameter circle
|
|---|
| 1909 | * \param CircleRadius radius of the parameter circle
|
|---|
| 1910 | * \param NewSphereCenter new center of a circumcircle
|
|---|
| 1911 | * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
|
|---|
| 1912 | * \param NormalVector normal vector
|
|---|
| 1913 | * \param SearchDirection search direction to make angle unique on return.
|
|---|
| 1914 | * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
|
|---|
| 1915 | */
|
|---|
| 1916 | double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
|
|---|
| 1917 | {
|
|---|
| 1918 | Vector helper;
|
|---|
| 1919 | double radius, alpha;
|
|---|
| 1920 |
|
|---|
| 1921 | helper.CopyVector(&NewSphereCenter);
|
|---|
| 1922 | // test whether new center is on the parameter circle's plane
|
|---|
| 1923 | if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
|
|---|
| 1924 | cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
|
|---|
| 1925 | helper.ProjectOntoPlane(&CirclePlaneNormal);
|
|---|
| 1926 | }
|
|---|
| 1927 | radius = helper.ScalarProduct(&helper);
|
|---|
| 1928 | // test whether the new center vector has length of CircleRadius
|
|---|
| 1929 | if (fabs(radius - CircleRadius) > HULLEPSILON)
|
|---|
| 1930 | cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
|
|---|
| 1931 | alpha = helper.Angle(&OldSphereCenter);
|
|---|
| 1932 | // make the angle unique by checking the halfplanes/search direction
|
|---|
| 1933 | if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
|
|---|
| 1934 | alpha = 2.*M_PI - alpha;
|
|---|
| 1935 | cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
|
|---|
| 1936 | radius = helper.Distance(&OldSphereCenter);
|
|---|
| 1937 | helper.ProjectOntoPlane(&NormalVector);
|
|---|
| 1938 | // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
|
|---|
| 1939 | if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
|
|---|
| 1940 | cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
|
|---|
| 1941 | return alpha;
|
|---|
| 1942 | } else {
|
|---|
| 1943 | cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
|
|---|
| 1944 | return 2.*M_PI;
|
|---|
| 1945 | }
|
|---|
| 1946 | };
|
|---|
| 1947 |
|
|---|
| 1948 |
|
|---|
| 1949 | /** This recursive function finds a third point, to form a triangle with two given ones.
|
|---|
| 1950 | * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
|
|---|
| 1951 | * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
|
|---|
| 1952 | * the center of the sphere is still fixed up to a single parameter. The band of possible values
|
|---|
| 1953 | * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
|
|---|
| 1954 | * us the "null" on this circle, the new center of the candidate point will be some way along this
|
|---|
| 1955 | * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
|
|---|
| 1956 | * by the normal vector of the base triangle that always points outwards by construction.
|
|---|
| 1957 | * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
|
|---|
| 1958 | * We construct the normal vector that defines the plane this circle lies in, it is just in the
|
|---|
| 1959 | * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
|
|---|
| 1960 | * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
|
|---|
| 1961 | * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
|
|---|
| 1962 | * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
|
|---|
| 1963 | * both.
|
|---|
| 1964 | * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
|
|---|
| 1965 | * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
|
|---|
| 1966 | * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
|
|---|
| 1967 | * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
|
|---|
| 1968 | * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
|
|---|
| 1969 | * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
|
|---|
| 1970 | * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points
|
|---|
| 1971 | * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line
|
|---|
| 1972 | * @param OptCandidate candidate reference on return
|
|---|
| 1973 | * @param OptCandidateCenter candidate's sphere center on return
|
|---|
| 1974 | * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
|
|---|
| 1975 | * @param RADIUS radius of sphere
|
|---|
| 1976 | * @param *LC LinkedCell structure with neighbouring atoms
|
|---|
| 1977 | */
|
|---|
| 1978 | // void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
|
|---|
| 1979 | // {
|
|---|
| 1980 | // atom *Walker = NULL;
|
|---|
| 1981 | // Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
|
|---|
| 1982 | // Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
|
|---|
| 1983 | // Vector OldSphereCenter; // center of the sphere defined by the three points of BaseTriangle
|
|---|
| 1984 | // Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
|
|---|
| 1985 | // Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
|
|---|
| 1986 | // Vector NewNormalVector; // normal vector of the Candidate's triangle
|
|---|
| 1987 | // Vector SearchDirection; // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
|
|---|
| 1988 | // Vector helper;
|
|---|
| 1989 | // LinkedAtoms *List = NULL;
|
|---|
| 1990 | // double CircleRadius; // radius of this circle
|
|---|
| 1991 | // double radius;
|
|---|
| 1992 | // double alpha, Otheralpha; // angles (i.e. parameter for the circle).
|
|---|
| 1993 | // double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
|
|---|
| 1994 | // int N[NDIM], Nlower[NDIM], Nupper[NDIM];
|
|---|
| 1995 | // atom *Candidate = NULL;
|
|---|
| 1996 | //
|
|---|
| 1997 | // cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
|
|---|
| 1998 | //
|
|---|
| 1999 | // cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
|
|---|
| 2000 | //
|
|---|
| 2001 | // // construct center of circle
|
|---|
| 2002 | // CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
|
|---|
| 2003 | // CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
|
|---|
| 2004 | // CircleCenter.Scale(0.5);
|
|---|
| 2005 | //
|
|---|
| 2006 | // // construct normal vector of circle
|
|---|
| 2007 | // CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
|
|---|
| 2008 | // CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
|
|---|
| 2009 | //
|
|---|
| 2010 | // // calculate squared radius of circle
|
|---|
| 2011 | // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
|
|---|
| 2012 | // if (radius/4. < RADIUS*RADIUS) {
|
|---|
| 2013 | // CircleRadius = RADIUS*RADIUS - radius/4.;
|
|---|
| 2014 | // CirclePlaneNormal.Normalize();
|
|---|
| 2015 | // cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
|
|---|
| 2016 | //
|
|---|
| 2017 | // // construct old center
|
|---|
| 2018 | // GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
|
|---|
| 2019 | // helper.CopyVector(&BaseTriangle->NormalVector); // normal vector ensures that this is correct center of the two possible ones
|
|---|
| 2020 | // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
|
|---|
| 2021 | // helper.Scale(sqrt(RADIUS*RADIUS - radius));
|
|---|
| 2022 | // OldSphereCenter.AddVector(&helper);
|
|---|
| 2023 | // OldSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2024 | // cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
|
|---|
| 2025 | //
|
|---|
| 2026 | // // test whether old center is on the band's plane
|
|---|
| 2027 | // if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
|
|---|
| 2028 | // cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
|
|---|
| 2029 | // OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
|
|---|
| 2030 | // }
|
|---|
| 2031 | // radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
|
|---|
| 2032 | // if (fabs(radius - CircleRadius) < HULLEPSILON) {
|
|---|
| 2033 | //
|
|---|
| 2034 | // // construct SearchDirection
|
|---|
| 2035 | // SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
|
|---|
| 2036 | // helper.CopyVector(&BaseLine->endpoints[0]->node->x);
|
|---|
| 2037 | // for(int i=0;i<3;i++) // just take next different endpoint
|
|---|
| 2038 | // if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
|
|---|
| 2039 | // helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
|
|---|
| 2040 | // }
|
|---|
| 2041 | // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards!
|
|---|
| 2042 | // SearchDirection.Scale(-1.);
|
|---|
| 2043 | // SearchDirection.ProjectOntoPlane(&OldSphereCenter);
|
|---|
| 2044 | // SearchDirection.Normalize();
|
|---|
| 2045 | // cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
|
|---|
| 2046 | // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
|
|---|
| 2047 | // cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
|
|---|
| 2048 | // }
|
|---|
| 2049 | //
|
|---|
| 2050 | // if (LC->SetIndexToVector(&CircleCenter)) { // get cell for the starting atom
|
|---|
| 2051 | // for(int i=0;i<NDIM;i++) // store indices of this cell
|
|---|
| 2052 | // N[i] = LC->n[i];
|
|---|
| 2053 | // cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2054 | // } else {
|
|---|
| 2055 | // cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
|
|---|
| 2056 | // return;
|
|---|
| 2057 | // }
|
|---|
| 2058 | // // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
|
|---|
| 2059 | // cout << Verbose(2) << "LC Intervals:";
|
|---|
| 2060 | // for (int i=0;i<NDIM;i++) {
|
|---|
| 2061 | // Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
|
|---|
| 2062 | // Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
|
|---|
| 2063 | // cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
|
|---|
| 2064 | // }
|
|---|
| 2065 | // cout << endl;
|
|---|
| 2066 | // for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
|
|---|
| 2067 | // for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
|
|---|
| 2068 | // for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
|
|---|
| 2069 | // List = LC->GetCurrentCell();
|
|---|
| 2070 | // cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2071 | // if (List != NULL) {
|
|---|
| 2072 | // for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
|---|
| 2073 | // Candidate = (*Runner);
|
|---|
| 2074 | //
|
|---|
| 2075 | // // check for three unique points
|
|---|
| 2076 | // if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
|
|---|
| 2077 | // cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
|
|---|
| 2078 | //
|
|---|
| 2079 | // // construct both new centers
|
|---|
| 2080 | // GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
|
|---|
| 2081 | // OtherNewSphereCenter.CopyVector(&NewSphereCenter);
|
|---|
| 2082 | //
|
|---|
| 2083 | // if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
|
|---|
| 2084 | // helper.CopyVector(&NewNormalVector);
|
|---|
| 2085 | // cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
|
|---|
| 2086 | // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
|
|---|
| 2087 | // if (radius < RADIUS*RADIUS) {
|
|---|
| 2088 | // helper.Scale(sqrt(RADIUS*RADIUS - radius));
|
|---|
| 2089 | // cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
|
|---|
| 2090 | // NewSphereCenter.AddVector(&helper);
|
|---|
| 2091 | // NewSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2092 | // cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
|
|---|
| 2093 | //
|
|---|
| 2094 | // helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
|
|---|
| 2095 | // OtherNewSphereCenter.AddVector(&helper);
|
|---|
| 2096 | // OtherNewSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2097 | // cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
|
|---|
| 2098 | //
|
|---|
| 2099 | // // check both possible centers
|
|---|
| 2100 | // alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
|
|---|
| 2101 | // Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
|
|---|
| 2102 | // alpha = min(alpha, Otheralpha);
|
|---|
| 2103 | // if (*ShortestAngle > alpha) {
|
|---|
| 2104 | // OptCandidate = Candidate;
|
|---|
| 2105 | // *ShortestAngle = alpha;
|
|---|
| 2106 | // if (alpha != Otheralpha)
|
|---|
| 2107 | // OptCandidateCenter->CopyVector(&NewSphereCenter);
|
|---|
| 2108 | // else
|
|---|
| 2109 | // OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
|
|---|
| 2110 | // cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
|
|---|
| 2111 | // } else {
|
|---|
| 2112 | // if (OptCandidate != NULL)
|
|---|
| 2113 | // cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
|
|---|
| 2114 | // else
|
|---|
| 2115 | // cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
|
|---|
| 2116 | // }
|
|---|
| 2117 | //
|
|---|
| 2118 | // } else {
|
|---|
| 2119 | // cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
|
|---|
| 2120 | // }
|
|---|
| 2121 | // } else {
|
|---|
| 2122 | // cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
|
|---|
| 2123 | // }
|
|---|
| 2124 | // } else {
|
|---|
| 2125 | // cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
|
|---|
| 2126 | // }
|
|---|
| 2127 | // }
|
|---|
| 2128 | // }
|
|---|
| 2129 | // }
|
|---|
| 2130 | // } else {
|
|---|
| 2131 | // cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
|
|---|
| 2132 | // }
|
|---|
| 2133 | // } else {
|
|---|
| 2134 | // cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
|
|---|
| 2135 | // }
|
|---|
| 2136 | //
|
|---|
| 2137 | // cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
|
|---|
| 2138 | // };
|
|---|
| 2139 |
|
|---|
| 2140 |
|
|---|
| 2141 | /** Checks whether the triangle consisting of the three atoms is already present.
|
|---|
| 2142 | * Searches for the points in Tesselation::PointsOnBoundary and checks their
|
|---|
| 2143 | * lines. If any of the three edges already has two triangles attached, false is
|
|---|
| 2144 | * returned.
|
|---|
| 2145 | * \param *out output stream for debugging
|
|---|
| 2146 | * \param *Candidates endpoints of the triangle candidate
|
|---|
| 2147 | * \return false - triangle invalid due to edge criteria, true - triangle may be added.
|
|---|
| 2148 | */
|
|---|
| 2149 | bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
|
|---|
| 2150 | LineMap::iterator FindLine;
|
|---|
| 2151 | PointMap::iterator FindPoint;
|
|---|
| 2152 |
|
|---|
| 2153 | *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
|
|---|
| 2154 | for (int i=0;i<3;i++) { // check through all endpoints
|
|---|
| 2155 | FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
|
|---|
| 2156 | if (FindPoint != PointsOnBoundary.end())
|
|---|
| 2157 | TPS[i] = FindPoint->second;
|
|---|
| 2158 | else
|
|---|
| 2159 | TPS[i] = NULL;
|
|---|
| 2160 | }
|
|---|
| 2161 |
|
|---|
| 2162 | // check lines
|
|---|
| 2163 | for (int i=0;i<3;i++)
|
|---|
| 2164 | if (TPS[i] != NULL)
|
|---|
| 2165 | for (int j=i;j<3;j++)
|
|---|
| 2166 | if (TPS[j] != NULL) {
|
|---|
| 2167 | FindLine = TPS[i]->lines.find(TPS[j]->node->nr);
|
|---|
| 2168 | if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) {
|
|---|
| 2169 | *out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->second->TrianglesCount << " triangles attached." << endl;
|
|---|
| 2170 | *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
|
|---|
| 2171 | return false;
|
|---|
| 2172 | }
|
|---|
| 2173 | }
|
|---|
| 2174 | *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
|
|---|
| 2175 | return true;
|
|---|
| 2176 | };
|
|---|
| 2177 |
|
|---|
| 2178 | /** This recursive function finds a third point, to form a triangle with two given ones.
|
|---|
| 2179 | * Note that this function is for the starting triangle.
|
|---|
| 2180 | * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
|
|---|
| 2181 | * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
|
|---|
| 2182 | * the center of the sphere is still fixed up to a single parameter. The band of possible values
|
|---|
| 2183 | * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
|
|---|
| 2184 | * us the "null" on this circle, the new center of the candidate point will be some way along this
|
|---|
| 2185 | * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
|
|---|
| 2186 | * by the normal vector of the base triangle that always points outwards by construction.
|
|---|
| 2187 | * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
|
|---|
| 2188 | * We construct the normal vector that defines the plane this circle lies in, it is just in the
|
|---|
| 2189 | * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
|
|---|
| 2190 | * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
|
|---|
| 2191 | * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
|
|---|
| 2192 | * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
|
|---|
| 2193 | * both.
|
|---|
| 2194 | * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
|
|---|
| 2195 | * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
|
|---|
| 2196 | * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
|
|---|
| 2197 | * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
|
|---|
| 2198 | * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
|
|---|
| 2199 | * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
|
|---|
| 2200 | * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
|
|---|
| 2201 | * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
|
|---|
| 2202 | * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
|
|---|
| 2203 | * @param BaseLine BoundaryLineSet with the current base line
|
|---|
| 2204 | * @param ThirdNode third atom to avoid in search
|
|---|
| 2205 | * @param OptCandidate candidate reference on return
|
|---|
| 2206 | * @param OptCandidateCenter candidate's sphere center on return
|
|---|
| 2207 | * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
|
|---|
| 2208 | * @param RADIUS radius of sphere
|
|---|
| 2209 | * @param *LC LinkedCell structure with neighbouring atoms
|
|---|
| 2210 | */
|
|---|
| 2211 | 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)
|
|---|
| 2212 | {
|
|---|
| 2213 | Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
|
|---|
| 2214 | Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
|
|---|
| 2215 | Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
|
|---|
| 2216 | Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
|
|---|
| 2217 | Vector NewNormalVector; // normal vector of the Candidate's triangle
|
|---|
| 2218 | Vector helper;
|
|---|
| 2219 | LinkedAtoms *List = NULL;
|
|---|
| 2220 | double CircleRadius; // radius of this circle
|
|---|
| 2221 | double radius;
|
|---|
| 2222 | double alpha, Otheralpha; // angles (i.e. parameter for the circle).
|
|---|
| 2223 | int N[NDIM], Nlower[NDIM], Nupper[NDIM];
|
|---|
| 2224 | atom *Candidate = NULL;
|
|---|
| 2225 |
|
|---|
| 2226 | cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
|
|---|
| 2227 |
|
|---|
| 2228 | cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
|
|---|
| 2229 |
|
|---|
| 2230 | // construct center of circle
|
|---|
| 2231 | CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
|
|---|
| 2232 | CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
|
|---|
| 2233 | CircleCenter.Scale(0.5);
|
|---|
| 2234 |
|
|---|
| 2235 | // construct normal vector of circle
|
|---|
| 2236 | CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
|
|---|
| 2237 | CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
|
|---|
| 2238 |
|
|---|
| 2239 | // calculate squared radius of circle
|
|---|
| 2240 | radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
|
|---|
| 2241 | if (radius/4. < RADIUS*RADIUS) {
|
|---|
| 2242 | CircleRadius = RADIUS*RADIUS - radius/4.;
|
|---|
| 2243 | CirclePlaneNormal.Normalize();
|
|---|
| 2244 | cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
|
|---|
| 2245 |
|
|---|
| 2246 | // test whether old center is on the band's plane
|
|---|
| 2247 | if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
|
|---|
| 2248 | cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
|
|---|
| 2249 | OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
|
|---|
| 2250 | }
|
|---|
| 2251 | radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
|
|---|
| 2252 | if (fabs(radius - CircleRadius) < HULLEPSILON) {
|
|---|
| 2253 |
|
|---|
| 2254 | // check SearchDirection
|
|---|
| 2255 | cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
|
|---|
| 2256 | if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
|
|---|
| 2257 | cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
|
|---|
| 2258 | }
|
|---|
| 2259 | // get cell for the starting atom
|
|---|
| 2260 | if (LC->SetIndexToVector(&CircleCenter)) {
|
|---|
| 2261 | for(int i=0;i<NDIM;i++) // store indices of this cell
|
|---|
| 2262 | N[i] = LC->n[i];
|
|---|
| 2263 | cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2264 | } else {
|
|---|
| 2265 | cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
|
|---|
| 2266 | return;
|
|---|
| 2267 | }
|
|---|
| 2268 | // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
|
|---|
| 2269 | cout << Verbose(2) << "LC Intervals:";
|
|---|
| 2270 | for (int i=0;i<NDIM;i++) {
|
|---|
| 2271 | Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
|
|---|
| 2272 | Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
|
|---|
| 2273 | cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
|
|---|
| 2274 | }
|
|---|
| 2275 | cout << endl;
|
|---|
| 2276 | for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
|
|---|
| 2277 | for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
|
|---|
| 2278 | for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
|
|---|
| 2279 | List = LC->GetCurrentCell();
|
|---|
| 2280 | //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2281 | if (List != NULL) {
|
|---|
| 2282 | for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
|---|
| 2283 | Candidate = (*Runner);
|
|---|
| 2284 |
|
|---|
| 2285 | // check for three unique points
|
|---|
| 2286 | if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
|
|---|
| 2287 | cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
|
|---|
| 2288 |
|
|---|
| 2289 | // construct both new centers
|
|---|
| 2290 | GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
|
|---|
| 2291 | OtherNewSphereCenter.CopyVector(&NewSphereCenter);
|
|---|
| 2292 |
|
|---|
| 2293 | if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
|
|---|
| 2294 | helper.CopyVector(&NewNormalVector);
|
|---|
| 2295 | cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
|
|---|
| 2296 | radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
|
|---|
| 2297 | if (radius < RADIUS*RADIUS) {
|
|---|
| 2298 | helper.Scale(sqrt(RADIUS*RADIUS - radius));
|
|---|
| 2299 | cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
|
|---|
| 2300 | NewSphereCenter.AddVector(&helper);
|
|---|
| 2301 | NewSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2302 | cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
|
|---|
| 2303 |
|
|---|
| 2304 | helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
|
|---|
| 2305 | OtherNewSphereCenter.AddVector(&helper);
|
|---|
| 2306 | OtherNewSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2307 | cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
|
|---|
| 2308 |
|
|---|
| 2309 | // check both possible centers
|
|---|
| 2310 | alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
|
|---|
| 2311 | Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
|
|---|
| 2312 | alpha = min(alpha, Otheralpha);
|
|---|
| 2313 | if (*ShortestAngle > alpha) {
|
|---|
| 2314 | OptCandidate = Candidate;
|
|---|
| 2315 | *ShortestAngle = alpha;
|
|---|
| 2316 | if (alpha != Otheralpha)
|
|---|
| 2317 | OptCandidateCenter->CopyVector(&NewSphereCenter);
|
|---|
| 2318 | else
|
|---|
| 2319 | OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
|
|---|
| 2320 | cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
|
|---|
| 2321 | } else {
|
|---|
| 2322 | if (OptCandidate != NULL)
|
|---|
| 2323 | cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
|
|---|
| 2324 | else
|
|---|
| 2325 | cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
|
|---|
| 2326 | }
|
|---|
| 2327 |
|
|---|
| 2328 | } else {
|
|---|
| 2329 | cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
|
|---|
| 2330 | }
|
|---|
| 2331 | } else {
|
|---|
| 2332 | cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
|
|---|
| 2333 | }
|
|---|
| 2334 | } else {
|
|---|
| 2335 | if (ThirdNode != NULL)
|
|---|
| 2336 | cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
|
|---|
| 2337 | else
|
|---|
| 2338 | cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
|
|---|
| 2339 | }
|
|---|
| 2340 | }
|
|---|
| 2341 | }
|
|---|
| 2342 | }
|
|---|
| 2343 | } else {
|
|---|
| 2344 | cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
|
|---|
| 2345 | }
|
|---|
| 2346 | } else {
|
|---|
| 2347 | if (ThirdNode != NULL)
|
|---|
| 2348 | cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
|
|---|
| 2349 | else
|
|---|
| 2350 | cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
|
|---|
| 2351 | }
|
|---|
| 2352 |
|
|---|
| 2353 | cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
|
|---|
| 2354 | };
|
|---|
| 2355 |
|
|---|
| 2356 | /** Finds the second point of starting triangle.
|
|---|
| 2357 | * \param *a first atom
|
|---|
| 2358 | * \param *Candidate pointer to candidate atom on return
|
|---|
| 2359 | * \param Oben vector indicating the outside
|
|---|
| 2360 | * \param Opt_Candidate reference to recommended candidate on return
|
|---|
| 2361 | * \param Storage[3] array storing angles and other candidate information
|
|---|
| 2362 | * \param RADIUS radius of virtual sphere
|
|---|
| 2363 | * \param *LC LinkedCell structure with neighbouring atoms
|
|---|
| 2364 | */
|
|---|
| 2365 | void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
|
|---|
| 2366 | {
|
|---|
| 2367 | cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
|
|---|
| 2368 | Vector AngleCheck;
|
|---|
| 2369 | double norm = -1., angle;
|
|---|
| 2370 | LinkedAtoms *List = NULL;
|
|---|
| 2371 | int N[NDIM], Nlower[NDIM], Nupper[NDIM];
|
|---|
| 2372 |
|
|---|
| 2373 | if (LC->SetIndexToAtom(a)) { // get cell for the starting atom
|
|---|
| 2374 | for(int i=0;i<NDIM;i++) // store indices of this cell
|
|---|
| 2375 | N[i] = LC->n[i];
|
|---|
| 2376 | } else {
|
|---|
| 2377 | cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
|
|---|
| 2378 | return;
|
|---|
| 2379 | }
|
|---|
| 2380 | // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
|
|---|
| 2381 | cout << Verbose(2) << "LC Intervals:";
|
|---|
| 2382 | for (int i=0;i<NDIM;i++) {
|
|---|
| 2383 | Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
|
|---|
| 2384 | Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
|
|---|
| 2385 | cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
|
|---|
| 2386 | }
|
|---|
| 2387 | cout << endl;
|
|---|
| 2388 | for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
|
|---|
| 2389 | for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
|
|---|
| 2390 | for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
|
|---|
| 2391 | List = LC->GetCurrentCell();
|
|---|
| 2392 | //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2393 | if (List != NULL) {
|
|---|
| 2394 | for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
|---|
| 2395 | Candidate = (*Runner);
|
|---|
| 2396 | // check if we only have one unique point yet ...
|
|---|
| 2397 | if (a != Candidate) {
|
|---|
| 2398 | cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
|
|---|
| 2399 | AngleCheck.CopyVector(&(Candidate->x));
|
|---|
| 2400 | AngleCheck.SubtractVector(&(a->x));
|
|---|
| 2401 | norm = AngleCheck.Norm();
|
|---|
| 2402 | // second point shall have smallest angle with respect to Oben vector
|
|---|
| 2403 | if (norm < RADIUS) {
|
|---|
| 2404 | angle = AngleCheck.Angle(&Oben);
|
|---|
| 2405 | if (angle < Storage[0]) {
|
|---|
| 2406 | //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
|
|---|
| 2407 | cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
|
|---|
| 2408 | Opt_Candidate = Candidate;
|
|---|
| 2409 | Storage[0] = AngleCheck.Angle(&Oben);
|
|---|
| 2410 | //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
|
|---|
| 2411 | } else {
|
|---|
| 2412 | cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
|
|---|
| 2413 | }
|
|---|
| 2414 | } else {
|
|---|
| 2415 | cout << "Refused due to Radius " << norm << endl;
|
|---|
| 2416 | }
|
|---|
| 2417 | }
|
|---|
| 2418 | }
|
|---|
| 2419 | }
|
|---|
| 2420 | }
|
|---|
| 2421 | cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
|
|---|
| 2422 | };
|
|---|
| 2423 |
|
|---|
| 2424 | /** Finds the starting triangle for find_non_convex_border().
|
|---|
| 2425 | * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
|
|---|
| 2426 | * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
|
|---|
| 2427 | * point are called.
|
|---|
| 2428 | * \param RADIUS radius of virtual rolling sphere
|
|---|
| 2429 | * \param *LC LinkedCell structure with neighbouring atoms
|
|---|
| 2430 | */
|
|---|
| 2431 | void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
|
|---|
| 2432 | {
|
|---|
| 2433 | cout << Verbose(1) << "Begin of Find_starting_triangle\n";
|
|---|
| 2434 | int i = 0;
|
|---|
| 2435 | LinkedAtoms *List = NULL;
|
|---|
| 2436 | atom* FirstPoint;
|
|---|
| 2437 | atom* SecondPoint;
|
|---|
| 2438 | atom* MaxAtom[NDIM];
|
|---|
| 2439 | double max_coordinate[NDIM];
|
|---|
| 2440 | Vector Oben;
|
|---|
| 2441 | Vector helper;
|
|---|
| 2442 | Vector Chord;
|
|---|
| 2443 | Vector SearchDirection;
|
|---|
| 2444 | Vector OptCandidateCenter;
|
|---|
| 2445 |
|
|---|
| 2446 | Oben.Zero();
|
|---|
| 2447 |
|
|---|
| 2448 | for (i = 0; i < 3; i++) {
|
|---|
| 2449 | MaxAtom[i] = NULL;
|
|---|
| 2450 | max_coordinate[i] = -1;
|
|---|
| 2451 | }
|
|---|
| 2452 |
|
|---|
| 2453 | // 1. searching topmost atom with respect to each axis
|
|---|
| 2454 | for (int i=0;i<NDIM;i++) { // each axis
|
|---|
| 2455 | LC->n[i] = LC->N[i]-1; // current axis is topmost cell
|
|---|
| 2456 | for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
|
|---|
| 2457 | for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
|
|---|
| 2458 | List = LC->GetCurrentCell();
|
|---|
| 2459 | //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 2460 | if (List != NULL) {
|
|---|
| 2461 | for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
|
|---|
| 2462 | cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
|
|---|
| 2463 | if ((*Runner)->x.x[i] > max_coordinate[i]) {
|
|---|
| 2464 | max_coordinate[i] = (*Runner)->x.x[i];
|
|---|
| 2465 | MaxAtom[i] = (*Runner);
|
|---|
| 2466 | }
|
|---|
| 2467 | }
|
|---|
| 2468 | } else {
|
|---|
| 2469 | cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
|
|---|
| 2470 | }
|
|---|
| 2471 | }
|
|---|
| 2472 | }
|
|---|
| 2473 |
|
|---|
| 2474 | cout << Verbose(2) << "Found maximum coordinates: ";
|
|---|
| 2475 | for (int i=0;i<NDIM;i++)
|
|---|
| 2476 | cout << i << ": " << *MaxAtom[i] << "\t";
|
|---|
| 2477 | cout << endl;
|
|---|
| 2478 | const int k = 1; // arbitrary choice
|
|---|
| 2479 | Oben.x[k] = 1.;
|
|---|
| 2480 | FirstPoint = MaxAtom[k];
|
|---|
| 2481 | cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
|
|---|
| 2482 |
|
|---|
| 2483 | // add first point
|
|---|
| 2484 | AddTrianglePoint(FirstPoint, 0);
|
|---|
| 2485 |
|
|---|
| 2486 | double ShortestAngle;
|
|---|
| 2487 | atom* Opt_Candidate = NULL;
|
|---|
| 2488 | 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.
|
|---|
| 2489 |
|
|---|
| 2490 | 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_...
|
|---|
| 2491 | SecondPoint = Opt_Candidate;
|
|---|
| 2492 | cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
|
|---|
| 2493 |
|
|---|
| 2494 | // add second point and first baseline
|
|---|
| 2495 | AddTrianglePoint(SecondPoint, 1);
|
|---|
| 2496 | AddTriangleLine(TPS[0], TPS[1], 0);
|
|---|
| 2497 |
|
|---|
| 2498 | helper.CopyVector(&(FirstPoint->x));
|
|---|
| 2499 | helper.SubtractVector(&(SecondPoint->x));
|
|---|
| 2500 | helper.Normalize();
|
|---|
| 2501 | Oben.ProjectOntoPlane(&helper);
|
|---|
| 2502 | Oben.Normalize();
|
|---|
| 2503 | helper.VectorProduct(&Oben);
|
|---|
| 2504 | ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
|
|---|
| 2505 |
|
|---|
| 2506 | Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
|
|---|
| 2507 | Chord.SubtractVector(&(SecondPoint->x));
|
|---|
| 2508 | double radius = Chord.ScalarProduct(&Chord);
|
|---|
| 2509 | double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
|
|---|
| 2510 | helper.CopyVector(&Oben);
|
|---|
| 2511 | helper.Scale(CircleRadius);
|
|---|
| 2512 | // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
|
|---|
| 2513 |
|
|---|
| 2514 | cout << Verbose(2) << "Looking for third point candidates \n";
|
|---|
| 2515 | // look in one direction of baseline for initial candidate
|
|---|
| 2516 | Opt_Candidate = NULL;
|
|---|
| 2517 | SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ...
|
|---|
| 2518 |
|
|---|
| 2519 | cout << Verbose(1) << "Looking for third point candidates ...\n";
|
|---|
| 2520 | Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
|
|---|
| 2521 | cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
|
|---|
| 2522 |
|
|---|
| 2523 | // add third point
|
|---|
| 2524 | AddTrianglePoint(Opt_Candidate, 2);
|
|---|
| 2525 |
|
|---|
| 2526 | // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
|
|---|
| 2527 |
|
|---|
| 2528 | // Finally, we only have to add the found further lines
|
|---|
| 2529 | AddTriangleLine(TPS[1], TPS[2], 1);
|
|---|
| 2530 | AddTriangleLine(TPS[0], TPS[2], 2);
|
|---|
| 2531 | // ... and triangles to the Maps of the Tesselation class
|
|---|
| 2532 | BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
|
|---|
| 2533 | AddTriangleToLines();
|
|---|
| 2534 | // ... and calculate its normal vector (with correct orientation)
|
|---|
| 2535 | OptCandidateCenter.Scale(-1.);
|
|---|
| 2536 | cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
|
|---|
| 2537 | BTS->GetNormalVector(OptCandidateCenter);
|
|---|
| 2538 | cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
|
|---|
| 2539 | cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
|
|---|
| 2540 | cout << Verbose(1) << "End of Find_starting_triangle\n";
|
|---|
| 2541 | };
|
|---|
| 2542 |
|
|---|
| 2543 | /** This function finds a triangle to a line, adjacent to an existing one.
|
|---|
| 2544 | * @param out output stream for debugging
|
|---|
| 2545 | * @param *mol molecule with Atom's and Bond's
|
|---|
| 2546 | * @param Line current baseline to search from
|
|---|
| 2547 | * @param T current triangle which \a Line is edge of
|
|---|
| 2548 | * @param RADIUS radius of the rolling ball
|
|---|
| 2549 | * @param N number of found triangles
|
|---|
| 2550 | * @param *filename filename base for intermediate envelopes
|
|---|
| 2551 | * @param *LC LinkedCell structure with neighbouring atoms
|
|---|
| 2552 | */
|
|---|
| 2553 | bool Tesselation::Find_next_suitable_triangle(ofstream *out,
|
|---|
| 2554 | molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
|
|---|
| 2555 | const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
|
|---|
| 2556 | {
|
|---|
| 2557 | cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
|
|---|
| 2558 | ofstream *tempstream = NULL;
|
|---|
| 2559 | char NumberName[255];
|
|---|
| 2560 |
|
|---|
| 2561 | atom* Opt_Candidate = NULL;
|
|---|
| 2562 | Vector OptCandidateCenter;
|
|---|
| 2563 |
|
|---|
| 2564 | Vector CircleCenter;
|
|---|
| 2565 | Vector CirclePlaneNormal;
|
|---|
| 2566 | Vector OldSphereCenter;
|
|---|
| 2567 | Vector SearchDirection;
|
|---|
| 2568 | Vector helper;
|
|---|
| 2569 | atom *ThirdNode = NULL;
|
|---|
| 2570 | double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
|
|---|
| 2571 | double radius, CircleRadius;
|
|---|
| 2572 |
|
|---|
| 2573 | cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
|
|---|
| 2574 | for (int i=0;i<3;i++)
|
|---|
| 2575 | if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
|
|---|
| 2576 | ThirdNode = T.endpoints[i]->node;
|
|---|
| 2577 |
|
|---|
| 2578 | // construct center of circle
|
|---|
| 2579 | CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
|
|---|
| 2580 | CircleCenter.AddVector(&Line.endpoints[1]->node->x);
|
|---|
| 2581 | CircleCenter.Scale(0.5);
|
|---|
| 2582 |
|
|---|
| 2583 | // construct normal vector of circle
|
|---|
| 2584 | CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
|
|---|
| 2585 | CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
|
|---|
| 2586 |
|
|---|
| 2587 | // calculate squared radius of circle
|
|---|
| 2588 | radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
|
|---|
| 2589 | if (radius/4. < RADIUS*RADIUS) {
|
|---|
| 2590 | CircleRadius = RADIUS*RADIUS - radius/4.;
|
|---|
| 2591 | CirclePlaneNormal.Normalize();
|
|---|
| 2592 | cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
|
|---|
| 2593 |
|
|---|
| 2594 | // construct old center
|
|---|
| 2595 | GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
|
|---|
| 2596 | helper.CopyVector(&T.NormalVector); // normal vector ensures that this is correct center of the two possible ones
|
|---|
| 2597 | radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
|
|---|
| 2598 | helper.Scale(sqrt(RADIUS*RADIUS - radius));
|
|---|
| 2599 | OldSphereCenter.AddVector(&helper);
|
|---|
| 2600 | OldSphereCenter.SubtractVector(&CircleCenter);
|
|---|
| 2601 | cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
|
|---|
| 2602 |
|
|---|
| 2603 | // construct SearchDirection
|
|---|
| 2604 | SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
|
|---|
| 2605 | helper.CopyVector(&Line.endpoints[0]->node->x);
|
|---|
| 2606 | helper.SubtractVector(&ThirdNode->x);
|
|---|
| 2607 | if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards!
|
|---|
| 2608 | SearchDirection.Scale(-1.);
|
|---|
| 2609 | SearchDirection.ProjectOntoPlane(&OldSphereCenter);
|
|---|
| 2610 | SearchDirection.Normalize();
|
|---|
| 2611 | cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
|
|---|
| 2612 | if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
|
|---|
| 2613 | cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
|
|---|
| 2614 | }
|
|---|
| 2615 |
|
|---|
| 2616 | // add third point
|
|---|
| 2617 | cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
|
|---|
| 2618 | Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
|
|---|
| 2619 |
|
|---|
| 2620 | } else {
|
|---|
| 2621 | cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
|
|---|
| 2622 | }
|
|---|
| 2623 |
|
|---|
| 2624 | if (Opt_Candidate == NULL) {
|
|---|
| 2625 | cerr << "WARNING: Could not find a suitable candidate." << endl;
|
|---|
| 2626 | return false;
|
|---|
| 2627 | }
|
|---|
| 2628 | cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl;
|
|---|
| 2629 |
|
|---|
| 2630 | // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
|
|---|
| 2631 | atom *AtomCandidates[3];
|
|---|
| 2632 | AtomCandidates[0] = Opt_Candidate;
|
|---|
| 2633 | AtomCandidates[1] = Line.endpoints[0]->node;
|
|---|
| 2634 | AtomCandidates[2] = Line.endpoints[1]->node;
|
|---|
| 2635 | bool flag = CheckPresenceOfTriangle(out, AtomCandidates);
|
|---|
| 2636 |
|
|---|
| 2637 | if (flag) { // if so, add
|
|---|
| 2638 | AddTrianglePoint(Opt_Candidate, 0);
|
|---|
| 2639 | AddTrianglePoint(Line.endpoints[0]->node, 1);
|
|---|
| 2640 | AddTrianglePoint(Line.endpoints[1]->node, 2);
|
|---|
| 2641 |
|
|---|
| 2642 | AddTriangleLine(TPS[0], TPS[1], 0);
|
|---|
| 2643 | AddTriangleLine(TPS[0], TPS[2], 1);
|
|---|
| 2644 | AddTriangleLine(TPS[1], TPS[2], 2);
|
|---|
| 2645 |
|
|---|
| 2646 | BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
|
|---|
| 2647 | AddTriangleToLines();
|
|---|
| 2648 |
|
|---|
| 2649 | OptCandidateCenter.Scale(-1.);
|
|---|
| 2650 | BTS->GetNormalVector(OptCandidateCenter);
|
|---|
| 2651 | OptCandidateCenter.Scale(-1.);
|
|---|
| 2652 |
|
|---|
| 2653 | cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
|
|---|
| 2654 | cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
|
|---|
| 2655 | } else { // else, yell and do nothing
|
|---|
| 2656 | cout << Verbose(1) << "This triangle consisting of ";
|
|---|
| 2657 | cout << *Opt_Candidate << ", ";
|
|---|
| 2658 | cout << *Line.endpoints[0]->node << " and ";
|
|---|
| 2659 | cout << *Line.endpoints[1]->node << " ";
|
|---|
| 2660 | cout << "is invalid!" << endl;
|
|---|
| 2661 | return false;
|
|---|
| 2662 | }
|
|---|
| 2663 |
|
|---|
| 2664 | if (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
|
|---|
| 2665 | sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
|
|---|
| 2666 | if (DoTecplotOutput) {
|
|---|
| 2667 | string NameofTempFile(tempbasename);
|
|---|
| 2668 | NameofTempFile.append(NumberName);
|
|---|
| 2669 | for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
|
|---|
| 2670 | NameofTempFile.erase(npos, 1);
|
|---|
| 2671 | NameofTempFile.append(TecplotSuffix);
|
|---|
| 2672 | cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
|
|---|
| 2673 | tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
|
|---|
| 2674 | write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
|
|---|
| 2675 | tempstream->close();
|
|---|
| 2676 | tempstream->flush();
|
|---|
| 2677 | delete(tempstream);
|
|---|
| 2678 | }
|
|---|
| 2679 |
|
|---|
| 2680 | if (DoRaster3DOutput) {
|
|---|
| 2681 | string NameofTempFile(tempbasename);
|
|---|
| 2682 | NameofTempFile.append(NumberName);
|
|---|
| 2683 | for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
|
|---|
| 2684 | NameofTempFile.erase(npos, 1);
|
|---|
| 2685 | NameofTempFile.append(Raster3DSuffix);
|
|---|
| 2686 | cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
|
|---|
| 2687 | tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
|
|---|
| 2688 | write_raster3d_file(out, tempstream, this, mol);
|
|---|
| 2689 | // include the current position of the virtual sphere in the temporary raster3d file
|
|---|
| 2690 | // make the circumsphere's center absolute again
|
|---|
| 2691 | helper.CopyVector(&Line.endpoints[0]->node->x);
|
|---|
| 2692 | helper.AddVector(&Line.endpoints[1]->node->x);
|
|---|
| 2693 | helper.Scale(0.5);
|
|---|
| 2694 | OptCandidateCenter.AddVector(&helper);
|
|---|
| 2695 | Vector *center = mol->DetermineCenterOfAll(out);
|
|---|
| 2696 | OptCandidateCenter.AddVector(center);
|
|---|
| 2697 | delete(center);
|
|---|
| 2698 | // and add to file plus translucency object
|
|---|
| 2699 | *tempstream << "# current virtual sphere\n";
|
|---|
| 2700 | *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n";
|
|---|
| 2701 | *tempstream << "2\n " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
|
|---|
| 2702 | *tempstream << "9\n terminating special property\n";
|
|---|
| 2703 | tempstream->close();
|
|---|
| 2704 | tempstream->flush();
|
|---|
| 2705 | delete(tempstream);
|
|---|
| 2706 | }
|
|---|
| 2707 | if (DoTecplotOutput || DoRaster3DOutput)
|
|---|
| 2708 | TriangleFilesWritten++;
|
|---|
| 2709 | }
|
|---|
| 2710 |
|
|---|
| 2711 | cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
|
|---|
| 2712 | return true;
|
|---|
| 2713 | };
|
|---|
| 2714 |
|
|---|
| 2715 | /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
|
|---|
| 2716 | * \param *out output stream for debugging
|
|---|
| 2717 | * \param *mol molecule structure with Atom's and Bond's
|
|---|
| 2718 | * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
|
|---|
| 2719 | * \param *LCList linked cell list of all atoms
|
|---|
| 2720 | * \param *filename filename prefix for output of vertex data
|
|---|
| 2721 | * \para RADIUS radius of the virtual sphere
|
|---|
| 2722 | */
|
|---|
| 2723 | void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
|
|---|
| 2724 | {
|
|---|
| 2725 | int N = 0;
|
|---|
| 2726 | bool freeTess = false;
|
|---|
| 2727 | *out << Verbose(1) << "Entering search for non convex hull. " << endl;
|
|---|
| 2728 | if (Tess == NULL) {
|
|---|
| 2729 | *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
|
|---|
| 2730 | Tess = new Tesselation;
|
|---|
| 2731 | freeTess = true;
|
|---|
| 2732 | }
|
|---|
| 2733 | bool freeLC = false;
|
|---|
| 2734 | LineMap::iterator baseline;
|
|---|
| 2735 | *out << Verbose(0) << "Begin of Find_non_convex_border\n";
|
|---|
| 2736 | bool flag = false; // marks whether we went once through all baselines without finding any without two triangles
|
|---|
| 2737 | bool failflag = false;
|
|---|
| 2738 |
|
|---|
| 2739 | if (LCList == NULL) {
|
|---|
| 2740 | LCList = new LinkedCell(mol, 2.*RADIUS);
|
|---|
| 2741 | freeLC = true;
|
|---|
| 2742 | }
|
|---|
| 2743 |
|
|---|
| 2744 | Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
|
|---|
| 2745 |
|
|---|
| 2746 | baseline = Tess->LinesOnBoundary.begin();
|
|---|
| 2747 | while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
|
|---|
| 2748 | if (baseline->second->TrianglesCount == 1) {
|
|---|
| 2749 | 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.
|
|---|
| 2750 | flag = flag || failflag;
|
|---|
| 2751 | if (!failflag)
|
|---|
| 2752 | cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
|
|---|
| 2753 | } else {
|
|---|
| 2754 | cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
|
|---|
| 2755 | }
|
|---|
| 2756 | N++;
|
|---|
| 2757 | baseline++;
|
|---|
| 2758 | if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
|
|---|
| 2759 | baseline = Tess->LinesOnBoundary.begin(); // restart if we reach end due to newly inserted lines
|
|---|
| 2760 | flag = false;
|
|---|
| 2761 | }
|
|---|
| 2762 | }
|
|---|
| 2763 | if (1) { //failflag) {
|
|---|
| 2764 | *out << Verbose(1) << "Writing final tecplot file\n";
|
|---|
| 2765 | if (DoTecplotOutput) {
|
|---|
| 2766 | string OutputName(filename);
|
|---|
| 2767 | OutputName.append(TecplotSuffix);
|
|---|
| 2768 | ofstream *tecplot = new ofstream(OutputName.c_str());
|
|---|
| 2769 | write_tecplot_file(out, tecplot, Tess, mol, -1);
|
|---|
| 2770 | tecplot->close();
|
|---|
| 2771 | delete(tecplot);
|
|---|
| 2772 | }
|
|---|
| 2773 | if (DoRaster3DOutput) {
|
|---|
| 2774 | string OutputName(filename);
|
|---|
| 2775 | OutputName.append(Raster3DSuffix);
|
|---|
| 2776 | ofstream *raster = new ofstream(OutputName.c_str());
|
|---|
| 2777 | write_raster3d_file(out, raster, Tess, mol);
|
|---|
| 2778 | raster->close();
|
|---|
| 2779 | delete(raster);
|
|---|
| 2780 | }
|
|---|
| 2781 | } else {
|
|---|
| 2782 | cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
|
|---|
| 2783 | }
|
|---|
| 2784 | if (freeTess)
|
|---|
| 2785 | delete(Tess);
|
|---|
| 2786 | if (freeLC)
|
|---|
| 2787 | delete(LCList);
|
|---|
| 2788 | *out << Verbose(0) << "End of Find_non_convex_border\n";
|
|---|
| 2789 | };
|
|---|
| 2790 |
|
|---|
| 2791 | /** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
|
|---|
| 2792 | * \param *out output stream for debugging
|
|---|
| 2793 | * \param *srcmol molecule to embed into
|
|---|
| 2794 | * \return *Vector new center of \a *srcmol for embedding relative to \a this
|
|---|
| 2795 | */
|
|---|
| 2796 | Vector* molecule::FindEmbeddingHole(ofstream *out, molecule *srcmol)
|
|---|
| 2797 | {
|
|---|
| 2798 | Vector *Center = new Vector;
|
|---|
| 2799 | Center->Zero();
|
|---|
| 2800 | // calculate volume/shape of \a *srcmol
|
|---|
| 2801 |
|
|---|
| 2802 | // find embedding holes
|
|---|
| 2803 |
|
|---|
| 2804 | // if more than one, let user choose
|
|---|
| 2805 |
|
|---|
| 2806 | // return embedding center
|
|---|
| 2807 | return Center;
|
|---|
| 2808 | };
|
|---|
| 2809 |
|
|---|
