Changeset 9d6308

Timestamp:

Jan 10, 2010, 7:24:02 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:

2793ef

Parents:

31a6af

Message:

Vector::GetIntersectionOfTwoLinesOnPlane() rewritten.

Use routine suggested by Eric Weisstein of Wolfram, which is basically what we did before the matrix attempt.

Signed-off-by: Frederik Heber <heber@tabletINS.(none)>

File:

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

src/vector.cpp

@@ -9 +9 @@
 #include "helpers.hpp"
 #include "info.hpp"
+#include "gslmatrix.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
@@ -293 +294 @@
 
 /** 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
@@ -307 +306 @@
 {
   Info FunctionInfo(__func__);
-  Vector Direction, OtherDirection;
-  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
-  gsl_vector *b = gsl_vector_alloc(NDIM);
-  gsl_vector *x = gsl_vector_alloc(NDIM);
-  gsl_permutation *perm = NULL;
-  int signum;
-
-  // construct both direction vectors
-  Direction.CopyVector(Line1b);
-  Direction.SubtractVector(Line1a);
-  if (Direction.IsZero())
+  Vector a;
+  Vector b;
+  Vector c;
+
+  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]);
+  }
+  if (fabs(M->Determinant()) < MYEPSILON) 
     return false;
-  OtherDirection.CopyVector(Line2b);
-  OtherDirection.SubtractVector(Line2a);
-  if (OtherDirection.IsZero())
-    return false;
-
-  // set vector
-  for (int i=0;i<NDIM;i++)
-    gsl_vector_set(b, i, Line1a->x[i]-Line2a->x[i]);
-  Log() << Verbose(1) << "b = " << endl;
-  gsl_vector_fprintf(stdout, b, "%g");
-
-  // set matrix
-  for (int i=0;i<NDIM;i++)
-    gsl_matrix_set(A, 0, i, -Direction.x[i]);
-  for (int i=0;i<NDIM;i++)
-    gsl_matrix_set(A, 1, i, OtherDirection.x[i]);
-  for (int i=0;i<NDIM;i++)
-    gsl_matrix_set(A, 2, i, 1.);
-  Log() << Verbose(1) << "A = " << endl;
-  gsl_matrix_fprintf(stdout, A, "%g");
-
-  // Solve A x=b for x
-  perm = gsl_permutation_alloc(NDIM);
-  gsl_linalg_LU_decomp(A, perm, &signum);
-  gsl_linalg_LU_solve(A, perm, b, x);
-  gsl_permutation_free(perm);
-  gsl_vector_free(b);
-  gsl_matrix_free(A);
-
-  Log() << Verbose(1) << "Solution is " << gsl_vector_get(x,0) << ", " << gsl_vector_get(x,1) << "." << endl;
-
-  CopyVector(&Direction);
-  Scale(gsl_vector_get(x,0));
+
+  // constuct a,b,c
+  a.CopyVector(Line1b);
+  a.SubtractVector(Line1a);
+  b.CopyVector(Line2b);
+  b.SubtractVector(Line2a);
+  c.CopyVector(Line2a);
+  c.SubtractVector(Line1a);
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1.CopyVector(&c);
+  temp1.VectorProduct(&b);
+  temp2.CopyVector(&a);
+  temp2.VectorProduct(&b);
+  s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+  Log() << Verbose(1) << "Factor s is " << s << "." << endl;
+
+  // construct intersection
+  CopyVector(&a);
+  Scale(s);
   AddVector(Line1a);
-  Log() << Verbose(1) << "INFO: First intersection is " << *this << "." << endl;
-
-  Vector test;
-  test.CopyVector(&OtherDirection);
-  test.Scale(gsl_vector_get(x,1));
-  test.AddVector(Line2a);
-  Log() << Verbose(1) << "INFO: Second intersection is " << test << "." << endl;
-  test.SubtractVector(this);
-
-  gsl_vector_free(x);
-
-  if (test.IsZero())
-    return true;
-  else
-    return false;
+  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
+
+  return true;
 };