Changeset a92c27 for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

May 25, 2016, 7:13:59 AM (9 years ago)

Author:

Frederik Heber <heber@…>

Children:

46e561

Parents:

4492f1

git-author:

Frederik Heber <heber@…> (07/12/14 16:19:46)

git-committer:

Frederik Heber <heber@…> (05/25/16 07:13:59)

Message:

Removed all IsZero() checking and special cases for the triangle idea.

File:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (3 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -697 +697 @@
       _referencepositions, _currentpositions);
 
-  if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
-    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
-    oldCenter.Normalize();
-    newCenter.Normalize();
-    if (!oldCenter.IsEqualTo(newCenter)) {
-      // calculate rotation axis and angle
-      Rotation.first = oldCenter;
-      Rotation.first.VectorProduct(newCenter);
-      Rotation.second = oldCenter.Angle(newCenter); // /(M_PI/2.);
-    } else {
-      // no rotation required anymore
-    }
+  ASSERT( !oldCenter.IsZero() && !newCenter.IsZero(),
+      "findPlaneAligningRotation() - either old "+toString(oldCenter)
+      +" or new center "+toString(newCenter)+" are zero.");
+  LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
+  if (!oldCenter.IsEqualTo(newCenter)) {
+    // calculate rotation axis and angle
+    Rotation.first = oldCenter;
+    Rotation.first.VectorProduct(newCenter);
+    Rotation.first.Normalize();
+    // construct reference vector to determine direction of rotation
+    const double sign = determineSignOfRotation(newCenter, oldCenter, Rotation.first);
+    Rotation.second = sign * oldCenter.Angle(newCenter);
   } else {
-    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
-    if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
-      // either oldCenter or newCenter (or both) is directly at origin
-      if (_currentpositions.indices.size() == 2) {
-        // line case
-        Vector oldPosition = _currentpositions.polygon[*_currentpositions.indices.begin()];
-        Vector newPosition = _referencepositions.polygon[0];
-        // check whether we need to rotate at all
-        if (!oldPosition.IsEqualTo(newPosition)) {
-          Rotation.first = oldPosition;
-          Rotation.first.VectorProduct(newPosition);
-          // orientation will fix the sign here eventually
-          Rotation.second = oldPosition.Angle(newPosition);
-        } else {
-          // no rotation required anymore
-        }
-      } else {
-        // triangle case
-        // both triangles/planes have same center, hence get axis by
-        // VectorProduct of Normals
-        Plane newplane(_referencepositions.polygon[0], _referencepositions.polygon[1], _referencepositions.polygon[2]);
-        VectorArray_t vectors;
-        for (IndexList_t::const_iterator iter = _currentpositions.indices.begin();
-            iter != _currentpositions.indices.end(); ++iter)
-          vectors.push_back(_currentpositions.polygon[*iter]);
-        Plane oldplane(vectors[0], vectors[1], vectors[2]);
-        Vector oldPosition = oldplane.getNormal();
-        Vector newPosition = newplane.getNormal();
-        // check whether we need to rotate at all
-        if (!oldPosition.IsEqualTo(newPosition)) {
-          Rotation.first = oldPosition;
-          Rotation.first.VectorProduct(newPosition);
-          Rotation.first.Normalize();
-
-          // construct reference vector to determine direction of rotation
-          const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
-          Rotation.second = sign * oldPosition.Angle(newPosition);
-          LOG(5, "DEBUG: Rotating plane normals by " << Rotation.second
-              << " around axis " << Rotation.first);
-        } else {
-          // else do nothing
-        }
-      }
-    } else {
-      // TODO: we can't do anything here, but this case needs to be dealt with when
-      // we have no ideal geometries anymore
-      if ((oldCenter-newCenter).Norm() > warn_amplitude)
-        ELOG(2, "oldCenter is " << oldCenter << ", yet newCenter is " << newCenter);
-      // else they are considered close enough
-      dontcheck = true;
-    }
+    // no rotation required anymore
   }
 
@@ -801 +751 @@
       << oldPosition.Norm() << ") and newPosition is " << newPosition << " length(: "
       << newPosition.Norm() << ")");
+
   if (!oldPosition.IsEqualTo(newPosition)) {
-    if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
-      // we rotate at oldCenter and around the radial direction, which is again given
-      // by oldCenter.
-      Rotation.first = oldCenter;
-      Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
-      LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
-          << Rotation.first << " as axis.");
-      oldPosition -= oldCenter;
-      newPosition -= oldCenter;
-      oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
-      newPosition = (newPosition - newPosition.Projection(Rotation.first));
-      LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
-    } else {
-      if (remainingnew.indices.size() == 2) {
-        // line situation
-        try {
-          Plane oldplane(oldPosition, oldCenter, newPosition);
-          Rotation.first = oldplane.getNormal();
-          LOG(6, "DEBUG: Plane is " << oldplane << " and normal is " << Rotation.first);
-        } catch (LinearDependenceException &e) {
-          LOG(6, "DEBUG: Vectors defining plane are linearly dependent.");
-          // oldPosition and newPosition are on a line, just flip when not equal
-          if (!oldPosition.IsEqualTo(newPosition)) {
-            Rotation.first.Zero();
-            Rotation.first.GetOneNormalVector(oldPosition);
-            LOG(6, "DEBUG: For flipping we use Rotation.first " << Rotation.first);
-            assert( Rotation.first.ScalarProduct(oldPosition) < std::numeric_limits<double>::epsilon()*1e4);
-  //              Rotation.second = M_PI;
-          } else {
-            LOG(6, "DEBUG: oldPosition and newPosition are equivalent.");
-          }
-        }
-      } else {
-        // triangle situation
-        Plane oldplane(remainingold.polygon[0], remainingold.polygon[1], remainingold.polygon[2]);
-        Rotation.first = oldplane.getNormal();
-        LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
-        oldPosition.ProjectOntoPlane(Rotation.first);
-        LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
-      }
-    }
+    // we rotate at oldCenter and around the radial direction, which is again given
+    // by oldCenter.
+    Rotation.first = oldCenter;
+    Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+    LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+        << Rotation.first << " as axis.");
+    oldPosition -= oldCenter;
+    newPosition -= oldCenter;
+    oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+    newPosition = (newPosition - newPosition.Projection(Rotation.first));
+    LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
+
     // construct reference vector to determine direction of rotation
     const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
@@ -915 +837 @@
 #ifndef NDEBUG
       // check: rotated "newCenter" should now equal oldCenter
-      {
+      // we don't check in case of two points as these lie on a great circle
+      // and the center cannot stably be recalculated. We may reactivate this
+      // when we calculate centers only once
+      if (oldSet.indices.size() > 2) {
         Vector oldCenter;
         Vector rotatednewCenter;