Changeset a1acc5 for molecuilder/src/vector.cpp

Timestamp:

Jan 8, 2010, 2:04:22 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

22b47e

Parents:

eda56a

Message:

Attempt to fix the tesselation::IsInnerPoint().

We try the IsInnerPoint() as follows:

1. Find nearest BoundaryPoints - working
2. Find Closest BoundaryLine's - working
3. Find closest Triangle that is well aligned (wrt to NormalVector and Distance) - unsure whether correctly working
4. Check whether alignment is on boundary or inside/outside - working
5. If on boundary, we check whether it's inside of triangle by intersecting with boundary lines - not working

Hence, we code a wrapper for GSL routines, to - finally - allow for solution of linear system of equations.

Signed-off-by: Frederik Heber <heber@…>

File:

• molecuilder/src/vector.cpp (modified) (8 diffs)

molecuilder/src/vector.cpp

@@ -8 +8 @@
 #include "defs.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
+#include "info.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>
 
 /************************************ Functions for class vector ************************************/
@@ -219 +225 @@
 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 +233 @@
   Direction.SubtractVector(Origin);
   Direction.Normalize();
-  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
   if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
@@ -236 +243 @@
   factor = helper.ScalarProduct(PlaneNormal)/factor;
   if (factor < MYEPSILON) { // Origin is in-plane
-    //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
+    Log() << Verbose(1) << "Origin of line is in-plane, simple." << endl;
     CopyVector(Origin);
     return true;
@@ -243 +250 @@
   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 +257 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
   } else {
@@ -299 +306 @@
 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;
+  Info FunctionInfo(__func__);
   Vector Direction, OtherDirection;
-  Vector AuxiliaryNormal;
-  Vector Distance;
-  const Vector *Normal = NULL;
-  Vector *ConstructedNormal = NULL;
-  bool FreeNormal = false;
+  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
-  Zero();
   Direction.CopyVector(Line1b);
   Direction.SubtractVector(Line1a);
@@ -318 +324 @@
     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;
+  // 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));
+  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 {
-    // 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;
+  else
+    return false;
 };