Changeset fcad4b

Timestamp:

Jan 11, 2010, 3:06:30 AM (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:

fee69b

Parents:

68a53b

Message:

InsideOutside unit test of tesselation is working correctly.

• FIX: BoundaryTriangleSet::GetIntersectionInsideTriangle() - don't need helper, just check whether CrossPoint is returned (and true) for all of the three sides.
• FIX: Tesselation::IsInnerPoint() - projection onto plane and stuff was nonsense, just take the Point ans Intersection which is on the plane anyway.
• FIX: Vector::GetIntersectionOfTwoLinesOnPlane() - coefficient MUST be zero (then vectors are coplanar), but parallelity check was missing. Also, we have to check whether s is in [0,1] in order to see whether we are inside the triangle side or outside.

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

Location:

src

Files:

• tesselation.cpp (modified) (4 diffs)
• vector.cpp (modified) (4 diffs)

src/tesselation.cpp

@@ -420 +420 @@
     if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
       CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-      helper.CopyVector(endpoints[(i+1)%3]->node->node);
-      helper.SubtractVector(endpoints[i%3]->node->node);
+      i++;
+    } else 
       break;
-    } else 
-      i++;
   } while (i<3);
   if (i==3) {
-    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
-  }
-  Log() << Verbose(1) << "INFO: Crosspoint is " << CrossPoint << "." << endl;
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
-  } else { // outside!
-    Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
   }
@@ -3241 +3234 @@
 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
@@ -3252 +3245 @@
     return false;
   } else {
-    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << "." << endl;
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
   }
 
@@ -3259 +3252 @@
   DistanceToCenter.CopyVector(&Center);
   DistanceToCenter.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
 
   // check whether we are on boundary
   if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
     // calculate whether inside of triangle
-    DistanceToCenter.ProjectOntoPlane(&result->NormalVector);
-    DistanceToCenter.AddVector(&Center);
-    Center.CopyVector(&DistanceToCenter);
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
     Center.SubtractVector(&result->NormalVector); // points towards MolCenter
     DistanceToCenter.AddVector(&result->NormalVector); // points outside

src/vector.cpp

@@ -319 +319 @@
     M->Set(3, i, Line2b->x[i]);
   }
-  if (fabs(M->Determinant()) < MYEPSILON) 
+  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;
+  }
+  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
+
 
   // constuct a,b,c
@@ -329 +339 @@
   c.CopyVector(Line2a);
   c.SubtractVector(Line1a);
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+
+  // check for parallelity
+  Vector parallel;
+  parallel.CopyVector(&a);
+  parallel.SubtractVector(&b);
+  if (parallel.NormSquared() < MYEPSILON) {
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    return false;
+  }
 
   // obtain s
@@ -337 +357 @@
   temp2.CopyVector(&a);
   temp2.VectorProduct(&b);
-  s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
-  Log() << Verbose(1) << "Factor s is " << s << "." << endl;
+  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
@@ -346 +370 @@
   Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
 
-  return true;
+  if ((s >=0 ) && (s<=1))
+    return true;
+  else
+    return false;
 };