Changeset 32d7e4

Timestamp:

Oct 19, 2014, 5:13:10 PM (10 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:

d45a5b

Parents:

d33bb8

git-author:

Frederik Heber <heber@…> (10/07/14 06:57:30)

git-committer:

Frederik Heber <heber@…> (10/19/14 17:13:10)

Message:

FIX: Tesselation::FlipBaseline() caclculated NormalVectors wrong.

• before we just took the average of the normal vector of both present triangles as guide. Now, we check the intersection of the fourth endpoint on the plane of the first triangle and decide by on which side of the three boundary lines the intersection lies whether to flip the normal vectors sign for the associated (new) triangle or not.

File:

• src/Tesselation/tesselation.cpp (modified) (7 diffs)

src/Tesselation/tesselation.cpp

@@ -1902 +1902 @@
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
+  Vector BaseLineNormal[2];
+  Vector OtherEndpoint[2]; // fourth point to either triangle
   int OldTriangleNrs[2], OldBaseLineNr;
   int i, m;
 
   // calculate NormalVector for later use
-  BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
     ELOG(1, "Less than two triangles are attached to this baseline!");
     return NULL;
   }
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-    LOG(5, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
-    BaseLineNormal += (runner->second->NormalVector);
-  }
-  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+  {
+    int i=0;
+    for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+      LOG(5, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+      OtherEndpoint[(i+1)%2] = runner->second->GetThirdEndpoint(Base)->node->getPosition();
+      BaseLineNormal[i++] = (runner->second->NormalVector);
+      ASSERT( i <= 2,
+          "Tesselation::FlipBaseline() - not exactly two triangles found.");
+    }
+  }
 
   // get the two triangles
@@ -1946 +1951 @@
 
   // index OldLines and OldPoints
+  // note that oldlines[0,1] belong to first triangle and hence first normal
+  // vector and oldlines[2,3] belong to second triangle and its normal vector
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
     for (int j = 0; j < 3; j++) // all of their endpoints and baselines
@@ -1955 +1962 @@
         OldPoints[m++] = runner->second->endpoints[j];
 
+  Vector BasePoints[2]; // endpoints of Base
+  if (OldLines[0]->ContainsBoundaryPoint(Base->endpoints[0])) {
+    BasePoints[0] = Base->endpoints[0]->node->getPosition();
+    BasePoints[1] = Base->endpoints[1]->node->getPosition();
+  } else {
+    BasePoints[0] = Base->endpoints[1]->node->getPosition();
+    BasePoints[1] = Base->endpoints[0]->node->getPosition();
+  }
   // check whether everything is in place to create new lines and triangles
   if (i < 4) {
@@ -1970 +1985 @@
       return NULL;
     }
+
+  // construct plane of first triangle for calculating normal vectors later
+  const Plane triangleplane = Base->triangles.begin()->second->getPlane();
+  // get fourth point projected down onto this plane
+  const Vector Intersection = triangleplane.getClosestPoint(OtherEndpoint[0]);
 
   // remove triangles and baseline removes itself
@@ -1994 +2014 @@
   LOG(3, "DEBUG: Created new baseline " << *NewLine << ".");
 
+  // Explanation for normal vector choice:
+  // Decisive for the normal vector of the new triangle is whether the fourth
+  // endpoint is with respect to the joining boundary line on one side or on
+  // the other with respect to the endpoint of the plane triangle that is not
+  // contained in the joining boundary line.
+
   // construct new triangles with flipped baseline
   i = -1;
@@ -2005 +2031 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
-    BTS->GetNormalVector(BaseLineNormal);
+    {
+      Line joiningline = makeLineThrough(
+          OldLines[0]->endpoints[0]->node->getPosition(),
+          OldLines[0]->endpoints[1]->node->getPosition());
+      // BasePoints[1] is not contained in OldLines[0], hence is third endpoint
+      // of plane triangle. BasePoints[0] is in the joining OldLines[0] and
+      // we check whether Intersection is on same side as BasePoints[1] or not.
+      const Vector pointinginside =
+          joiningline.getClosestPoint(BasePoints[1]) - BasePoints[1];
+      const Vector pointingtointersection =
+          joiningline.getClosestPoint(Intersection) - Intersection;
+      const double sign_of_intersection =
+          pointingtointersection.ScalarProduct(pointinginside);
+      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[0]
+          << " is " << sign_of_intersection << ".");
+      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
+      LOG(4, "DEBUG: Opposite normal direction is "
+          << sign_of_normal * BaseLineNormal[0] << ".");
+      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
+    }
     AddTesselationTriangle(OldTriangleNrs[0]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
@@ -2013 +2058 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
-    BTS->GetNormalVector(BaseLineNormal);
+    {
+      Line joiningline = makeLineThrough(
+          OldLines[1]->endpoints[0]->node->getPosition(),
+          OldLines[1]->endpoints[1]->node->getPosition());
+      // BasePoints[0] is not contained in OldLines[1], hence is third endpoint
+      // of plane triangle. BasePoints[1] is in th1e joining OldLines[1] and
+      // we check whether Intersection is on same side as BasePoints[0] or not.
+      const Vector pointinginside =
+          joiningline.getClosestPoint(BasePoints[0]) - BasePoints[0];
+      const Vector pointingtointersection =
+          joiningline.getClosestPoint(Intersection) - Intersection;
+      const double sign_of_intersection =
+          pointingtointersection.ScalarProduct(pointinginside);
+      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[1]
+          << " is " << sign_of_intersection << ".");
+      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
+      LOG(4, "DEBUG: Opposite normal direction is "
+          << sign_of_normal * BaseLineNormal[0] << ".");
+      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
+    }
     AddTesselationTriangle(OldTriangleNrs[1]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");