Changeset 46670d

Timestamp:

Aug 8, 2009, 7:19:23 PM (15 years ago)

Author:

Frederik Heber <heber@…>

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:

8bb475f

Parents:

055861

Message:

BUGFIXES to some vector functions.

• Vector::MakeNormalVector() ... Projection lacks ../Norm()/Norm() and hence produced wrong vector.
• Vector::GetIntersectionWithPlane() - checks for uniqueness and Bronstein formula removed and correct one implemented.
• Vector::GetIntersectionOfTwoLinesOnPlane() - Bronstein formula removed, correct procedure implemented.
• Vector::IsNull() is const function.

Location:

src

Files:

• vector.cpp (modified) (4 diffs)
• vector.hpp (modified) (2 diffs)

src/vector.cpp

@@ -213 +213 @@
  * \return true -  \a this contains intersection point on return, false - line is parallel to plane
  */
-bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *LineVector, Vector *LineVector2)
+bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *Origin, Vector *LineVector)
 {
   double factor;
-  Vector Direction;
+  Vector Direction, helper;
 
   // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
-  Direction.CopyVector(LineVector2);
-  Direction.SubtractVector(LineVector);
-  if (Direction.ScalarProduct(PlaneNormal) < MYEPSILON)
+  Direction.CopyVector(LineVector);
+  Direction.SubtractVector(Origin);
+  factor = Direction.ScalarProduct(PlaneNormal);
+  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
+    *out << Verbose(2) << "WARNING: Line is parallel to plane, no intersection." << endl;
     return false;
-  factor = LineVector->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
+  }
+  helper.CopyVector(PlaneOffset);
+  helper.SubtractVector(LineVector);
+  factor = helper.ScalarProduct(PlaneNormal)/factor;
+  //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
   Direction.Scale(factor);
   CopyVector(LineVector);
-  SubtractVector(&Direction);
+  AddVector(&Direction);
 
   // test whether resulting vector really is on plane
-  Direction.CopyVector(this);
-  Direction.SubtractVector(PlaneOffset);
-  if (Direction.ScalarProduct(PlaneNormal) > MYEPSILON)
+  helper.CopyVector(this);
+  helper.SubtractVector(PlaneOffset);
+  if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
+    *out << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
-  else
+  } else {
+    *out << Verbose(2) << "WARNING: Intersection point " << *this << " is not on plane." << endl;
     return false;
+  }
 };
 
@@ -244 +253 @@
  * \param *Line2a first vector of second line
  * \param *Line2b second vector of second line
+ * \param *PlaneNormal normal of plane, is supplemental/arbitrary
  * \return true - \a this will contain the intersection on return, false - lines are parallel
  */
-bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b)
-{
-  double k1,a1,h1,b1,k2,a2,h2,b2;
-  // equation for line 1
-  if (fabs(Line1a->x[0] - Line2a->x[0]) < MYEPSILON) {
-    k1 = 0;
-    h1 = 0;
-  } else {
-    k1 = (Line1a->x[1] - Line2a->x[1])/(Line1a->x[0] - Line2a->x[0]);
-    h1 = (Line1a->x[2] - Line2a->x[2])/(Line1a->x[0] - Line2a->x[0]);
-  }
-  a1 = 0.5*((Line1a->x[1] + Line2a->x[1]) - k1*(Line1a->x[0] + Line2a->x[0]));
-  b1 = 0.5*((Line1a->x[2] + Line2a->x[2]) - h1*(Line1a->x[0] + Line2a->x[0]));
-
-  // equation for line 2
-  if (fabs(Line2a->x[0] - Line2a->x[0]) < MYEPSILON) {
-    k1 = 0;
-    h1 = 0;
-  } else {
-    k1 = (Line2a->x[1] - Line2a->x[1])/(Line2a->x[0] - Line2a->x[0]);
-    h1 = (Line2a->x[2] - Line2a->x[2])/(Line2a->x[0] - Line2a->x[0]);
-  }
-  a1 = 0.5*((Line2a->x[1] + Line2a->x[1]) - k1*(Line2a->x[0] + Line2a->x[0]));
-  b1 = 0.5*((Line2a->x[2] + Line2a->x[2]) - h1*(Line2a->x[0] + Line2a->x[0]));
-
-  // calculate cross point
-  if (fabs((a1-a2)*(h1-h2) - (b1-b2)*(k1-k2)) < MYEPSILON) {
-    x[0] = (a2-a1)/(k1-k2);
-    x[1] = (k1*a2-k2*a1)/(k1-k2);
-    x[2] = (h1*b2-h2*b1)/(h1-h2);
-    *out << Verbose(4) << "Lines do intersect at " << this << "." << endl;
-    return true;
-  } else {
-    *out << Verbose(4) << "Lines do not intersect." << endl;
-    return false;
-  }
+bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *PlaneNormal)
+{
+  double factor1, factor2;
+  Vector helper, Line, LineNormal, *OtherNormal = NULL;
+  const Vector *Normal;
+  bool result = false;
+
+  // create Plane normal vector
+  if (PlaneNormal == NULL) {
+    OtherNormal = new Vector(0.,0.,0.);
+    if (!OtherNormal->MakeNormalVector(Line1a, Line1b, Line2a))
+      if (!OtherNormal->MakeNormalVector(Line1a, Line1b, Line2b)) {
+        *out << Verbose(1) << "ERROR: GetIntersectionOfTwoLinesOnPlane() cannot create a normal of the plane, everything is linear dependent." << endl;
+        return false;
+      }
+    Normal = OtherNormal;
+  } else
+    Normal = PlaneNormal;
+  *out << Verbose(3) << "INFO: Normal of plane is " << *Normal << "." << endl;
+
+  // create normal vector to one line
+  Line.CopyVector(Line1b);
+  Line.SubtractVector(Line1a);
+  LineNormal.MakeNormalVector(&Line, Normal);
+  *out << Verbose(3) << "INFO: Normal of first line is " << LineNormal << "." << endl;
+
+  // check if lines are parallel
+  helper.CopyVector(Line2b);
+  helper.SubtractVector(Line2a);
+  if (fabs(helper.ScalarProduct(&LineNormal)) < MYEPSILON) {
+    *out << Verbose(1) << "Lines " << helper << " and " << Line << " are parallel, no cross point!" << endl;
+    result = false;
+  } else {  
+    helper.CopyVector(Line2a);
+    helper.SubtractVector(Line1a);
+    factor1 = helper.ScalarProduct(&LineNormal);
+    helper.CopyVector(Line2b);
+    helper.SubtractVector(Line1a);
+    factor2 = helper.ScalarProduct(&LineNormal);
+    if (fabs(factor2) > MYEPSILON) {
+      CopyVector(Line2a);
+      helper.Scale(factor1/factor2);
+      AddVector(&helper);
+      result = true;
+    } else {
+      Zero();
+      result = false;
+    }
+  }
+
+  if (OtherNormal != NULL)
+    delete(OtherNormal);
+
+  return result;
 };
 
@@ -352 +380 @@
  * @return true - vector is zero, false - vector is not
  */
-bool Vector::IsNull()
+bool Vector::IsNull() const
 {
   return (fabs(x[0]+x[1]+x[2]) < MYEPSILON);
@@ -692 +720 @@
 {
   bool result = false;
+  double factor = y1->Projection(this)/y1->Norm()/y1->Norm();
   Vector x1;
   x1.CopyVector(y1);
-  x1.Scale(x1.Projection(this));
+  x1.Scale(factor);
   SubtractVector(&x1);
   for (int i=NDIM;i--;)

src/vector.hpp

@@ -28 +28 @@
   double NormSquared() const;
   double Angle(const Vector *y) const;
-  bool IsNull();
+  bool IsNull() const;
 
   void AddVector(const Vector *y);
@@ -51 +51 @@
   void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors);
   double CutsPlaneAt(Vector *A, Vector *B, Vector *C);
-  bool GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *LineVector, Vector *LineVector2);
-  bool GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b);
+  bool GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *Origin, Vector *LineVector);
+  bool GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *Normal = NULL);
   bool GetOneNormalVector(const Vector *x1);
   bool MakeNormalVector(const Vector *y1);