Changeset b8d1aeb for src/vector.cpp

Timestamp:

Feb 2, 2010, 11:38:06 AM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

7ba324

Parents:

9fe36b (diff), 2ededc2 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'MenuRefactoring' into QT4Refactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/builder.cpp
molecuilder/src/unittests/Makefile.am

File:

• src/vector.cpp (modified) (12 diffs)

src/vector.cpp

@@ -8 +8 @@
 #include "defs.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
+#include "info.hpp"
+#include "gslmatrix.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
+#include "memoryallocator.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
+
+#include <cassert>
 
 /************************************ Functions for class vector ************************************/
@@ -215 +224 @@
  * \param *Origin first vector of line
  * \param *LineVector second vector of line
- * \return true -  \a this contains intersection point on return, false - line is parallel to plane
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
  */
 bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
 {
+  Info FunctionInfo(__func__);
   double factor;
   Vector Direction, helper;
@@ -226 +236 @@
   Direction.SubtractVector(Origin);
   Direction.Normalize();
-  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
+  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
-  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
-    eLog() << Verbose(2) << "Line is parallel to plane, no intersection." << endl;
+  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
+    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
     return false;
   }
@@ -235 +246 @@
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (factor < MYEPSILON) { // Origin is in-plane
-    //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
+  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
+    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
     CopyVector(Origin);
     return true;
@@ -243 +254 @@
   Direction.Scale(factor);
   CopyVector(Origin);
-  //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 
@@ -250 +261 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
     return true;
   } else {
@@ -286 +297 @@
 
 /** Calculates the intersection of the two lines that are both on the same plane.
- * We construct auxiliary plane with its vector normal to one line direction and the PlaneNormal, then a vector
- * from the first line's offset onto the plane. Finally, scale by factor is 1/cos(angle(line1,line2..)) = 1/SP(...), and
- * project onto the first line's direction and add its offset.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html 
  * \param *out output stream for debugging
  * \param *Line1a first vector of first line
@@ -299 +308 @@
 bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
 {
-  bool result = true;
-  Vector Direction, OtherDirection;
-  Vector AuxiliaryNormal;
-  Vector Distance;
-  const Vector *Normal = NULL;
-  Vector *ConstructedNormal = NULL;
-  bool FreeNormal = false;
-
-  // construct both direction vectors
-  Zero();
-  Direction.CopyVector(Line1b);
-  Direction.SubtractVector(Line1a);
-  if (Direction.IsZero())
+  Info FunctionInfo(__func__);
+
+  GSLMatrix *M = new GSLMatrix(4,4);
+
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a->x[i]);
+    M->Set(1, i, Line1b->x[i]);
+    M->Set(2, i, Line2a->x[i]);
+    M->Set(3, i, Line2b->x[i]);
+  }
+  
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
     return false;
-  OtherDirection.CopyVector(Line2b);
-  OtherDirection.SubtractVector(Line2a);
-  if (OtherDirection.IsZero())
+  }
+  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
+
+
+  // constuct a,b,c
+  Vector a;
+  Vector b;
+  Vector c;
+  Vector d;
+  a.CopyVector(Line1b);
+  a.SubtractVector(Line1a);
+  b.CopyVector(Line2b);
+  b.SubtractVector(Line2a);
+  c.CopyVector(Line2a);
+  c.SubtractVector(Line1a);
+  d.CopyVector(Line2b);
+  d.SubtractVector(Line1b);
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   return false;
+  }
+
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel.CopyVector(Line1a);
+    parallel.SubtractVector(Line2a);
+    factor = parallel.ScalarProduct(&a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      CopyVector(Line2a);
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return true;
+    } else {
+      parallel.CopyVector(Line1a);
+      parallel.SubtractVector(Line2b);
+      factor = parallel.ScalarProduct(&a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        CopyVector(Line2b);
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return true;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    Zero();
     return false;
-
-  Direction.Normalize();
-  OtherDirection.Normalize();
-
-  //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
-
-  if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
-    if ((Line1a == Line2a) || (Line1a == Line2b))
-      CopyVector(Line1a);
-    else if ((Line1b == Line2b) || (Line1b == Line2b))
-        CopyVector(Line1b);
-    else
-      return false;
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
-    return true;
-  } else {
-    // check whether we have a plane normal vector
-    if (PlaneNormal == NULL) {
-      ConstructedNormal = new Vector;
-      ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
-      Normal = ConstructedNormal;
-      FreeNormal = true;
-    } else
-      Normal = PlaneNormal;
-
-    AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
-    //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
-
-    Distance.CopyVector(Line2a);
-    Distance.SubtractVector(Line1a);
-    //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
-    if (Distance.IsZero()) {
-      // offsets are equal, match found
-      CopyVector(Line1a);
-      result = true;
-    } else {
-      CopyVector(Distance.Projection(&AuxiliaryNormal));
-      //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
-      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
-      //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
-      Scale(1./(factor*factor));
-      //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
-      CopyVector(Projection(&Direction));
-      //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
-      if (this->IsZero())
-        result = false;
-      else
-        result = true;
-      AddVector(Line1a);
-    }
-
-    if (FreeNormal)
-      delete(ConstructedNormal);
-  }
-  if (result)
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
-
-  return result;
+  }
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1.CopyVector(&c);
+  temp1.VectorProduct(&b);
+  temp2.CopyVector(&a);
+  temp2.VectorProduct(&b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+
+  // construct intersection
+  CopyVector(&a);
+  Scale(s);
+  AddVector(Line1a);
+  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
+
+  return true;
 };
 
@@ -480 +505 @@
   else
     return false;
+};
+
+/** Checks whether vector is normal to \a *normal.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsEqualTo(const Vector * const a) const
+{
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a->x[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };
 
@@ -626 +664 @@
   return *x;
 };
+
+Vector& Vector::operator=(const Vector& src) {
+  CopyVector(src);
+  return *this;
+}
+
+double& Vector::operator[](int i){
+  assert(i<NDIM && "Invalid Vector dimension requested");
+  return x[i];
+}
 
 /** Prints a 3dim vector.
@@ -1040 +1088 @@
 void Vector::CopyVector(const Vector * const y)
 {
-  for (int i=NDIM;i--;)
-    this->x[i] = y->x[i];
+  // check for self assignment
+  if(y!=this){
+    for (int i=NDIM;i--;)
+      this->x[i] = y->x[i];
+  }
 }
 
@@ -1049 +1100 @@
 void Vector::CopyVector(const Vector &y)
 {
-  for (int i=NDIM;i--;)
-    this->x[i] = y.x[i];
+  // check for self assignment
+  if(&y!=this) {
+    for (int i=NDIM;i--;)
+      this->x[i] = y.x[i];
+  }
 }