Changeset 71910a for molecuilder/src/tesselation.cpp

Timestamp:

Apr 7, 2010, 3:45:38 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

Children:

0f55b2

Parents:

770138

Message:

Made data internal data-structure of vector class private

• Replaced occurences of access to internals with operator
• moved Vector-class into LinAlg-Module
• Reworked Vector to allow clean modularization
• Added Plane class to describe arbitrary planes in 3d space

File:

• molecuilder/src/tesselation.cpp (modified) (22 diffs)

molecuilder/src/tesselation.cpp

@@ -15 +15 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "vector_ops.hpp"
 #include "verbose.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 class molecule;
@@ -258 +261 @@
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
-      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
+      helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
@@ -402 +405 @@
         Info FunctionInfo(__func__);
   // get normal vector
-  NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
+  NormalVector = Plane(*(endpoints[0]->node->node),
+                       *(endpoints[1]->node->node),
+                       *(endpoints[2]->node->node)).getNormal();
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
@@ -428 +433 @@
   Vector helper;
 
-  if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
+  try {
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
+  }
+  catch (LinearDependenceException &excp) {
+    Log() << Verbose(1) << excp;
     eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
     return false;
@@ -450 +459 @@
   int i=0;
   do {
-    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
+    try {
+      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
+                                                    *(endpoints[(i+1)%3]->node->node),
+                                                    *(endpoints[(i+2)%3]->node->node),
+                                                    *Intersection);
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
@@ -462 +475 @@
       }
       i++;
-    } else 
+    } catch (LinearDependenceException &excp){
       break;
+    }
   } while (i<3);
   if (i==3) {
@@ -492 +506 @@
   Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
   GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
+  try {
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
+  }
+  catch (LinearDependenceException &excp) {
     ClosestPoint->CopyVector(x);
   }
@@ -518 +535 @@
     Direction.SubtractVector(endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
+
+    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+
     CrossDirection[i].CopyVector(&CrossPoint[i]);
     CrossDirection[i].SubtractVector(&InPlane);
@@ -731 +750 @@
   int counter=0;
   for (; Runner[2] != endpoints.end(); ) {
-    TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);
+    TemporaryNormal = Plane(*((*Runner[0])->node->node),
+                            *((*Runner[1])->node->node),
+                            *((*Runner[2])->node->node)).getNormal();
     for (int i=0;i<3;i++) // increase each of them
       Runner[i]++;
@@ -1220 +1241 @@
       // 2. next, we have to check whether all points reside on only one side of the triangle
       // 3. construct plane vector
-      PlaneVector.MakeNormalVector(A->second->node->node,
-          baseline->second.first->second->node->node,
-          baseline->second.second->second->node->node);
+      PlaneVector = Plane(*(A->second->node->node),
+                          *(baseline->second.first->second->node->node),
+                          *(baseline->second.second->second->node->node)).getNormal();
       Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
       // 4. loop over all points
@@ -1391 +1412 @@
         // vector in propagation direction (out of triangle)
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
+        PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
         TempVector.CopyVector(&CenterVector);
         TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
@@ -1448 +1469 @@
             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
             flag = true;
-            VirtualNormalVector.MakeNormalVector(baseline->second->endpoints[0]->node->node, baseline->second->endpoints[1]->node->node, target->second->node->node);
+            VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
+                                        *(baseline->second->endpoints[1]->node->node),
+                                        *(target->second->node->node)).getNormal();
             TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
             TempVector.AddVector(baseline->second->endpoints[1]->node->node);
@@ -2063 +2086 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if ((*Runner)->node->x[i] > maxCoordinate[i]) {
+            if ((*Runner)->node->at(i) > maxCoordinate[i]) {
               Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
-              maxCoordinate[i] = (*Runner)->node->x[i];
+              maxCoordinate[i] = (*Runner)->node->at(i);
               MaxPoint[i] = (*Runner);
             }
@@ -2083 +2106 @@
   for (int k=0;k<NDIM;k++) {
     NormalVector.Zero();
-    NormalVector.x[k] = 1.;
+    NormalVector[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
@@ -2117 +2140 @@
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -2336 +2359 @@
 
     // construct SearchDirection and an "outward pointer"
-    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
     helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
@@ -2979 +3002 @@
                   Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
 
-                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
-                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
-                  ) {
+                  try {
+                    NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
+                                              *(CandidateLine.BaseLine->endpoints[1]->node->node),
+                                              *(Candidate->node)).getNormal();
+
+                    if(!fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
+                      throw LinearDependenceException(__FILE__,__LINE__);
+
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                     radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
@@ -3043 +3071 @@
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
                     }
-                  } else {
+                  }
+                  catch (LinearDependenceException &excp){
+                    Log() << Verbose(1) << excp;
                     Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                   }
@@ -3565 +3595 @@
   Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
 
@@ -3633 +3663 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    helper = Plane(*((*TesselA)->node),
+                   *((*TesselB)->node),
+                   *((*TesselC)->node)).getNormal();
     Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
     counter++;
@@ -3670 +3702 @@
   Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
 
@@ -4002 +4034 @@
           OrthogonalVector.CopyVector((*MiddleNode)->node);
           OrthogonalVector.SubtractVector(&OldPoint);
-          OrthogonalVector.MakeNormalVector(&Reference);
+          OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
           //if (angle < M_PI)  // no wrong-sided triangles, please?