Changeset 7c6712 for src/boundary.cpp

Timestamp:

Dec 19, 2008, 4:21:11 PM (16 years ago)

Author:

Christian Neuen <neuen@…>

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:

196a5a, 44fd95

Parents:

e9fa06f

Message:

Switched choice process to Ball angle completely.
Problem remains.
Same position?

File:

• src/boundary.cpp (modified) (10 diffs)

src/boundary.cpp

@@ -1415 +1415 @@
 ;
 
+/**
+ * Function returns center of sphere with RADIUS, which rests on points a, b, c
+ * @param Center this vector will be used for return
+ * @param a vector first point of triangle
+ * @param b vector second point of triangle
+ * @param c vector third point of triangle
+ * @param Direction vector indicates up/down
+ * @param Halfplaneindicator double indicates whether Direction is up or down
+ * @param alpha double angle at a
+ * @param beta double, angle at b
+ * @param gamma, double, angle at c
+ * @param Radius, double
+ * @param Umkreisradius double radius of circumscribing circle
+ */
+
+  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, double HalfplaneIndicator,
+      double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+  {
+    Vector TempNormal, helper;
+    double Restradius;
+
+    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+    Center->Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    // Here we calculated center of circumscribing circle, using barycentric coordinates
+
+    TempNormal.CopyVector(&a);
+    TempNormal.SubtractVector(&b);
+    helper.CopyVector(&a);
+    helper.SubtractVector(&c);
+    TempNormal.VectorProduct(&helper);
+    if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+      {
+        TempNormal.Scale(-1);
+      }
+
+    TempNormal.Normalize();
+    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+    TempNormal.Scale(Restradius);
+
+    Center->AddVector(&TempNormal);
+  }
+  ;
+
+
 /** This recursive function finds a third point, to form a triangle with two given ones.
  * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
@@ -1429 +1473 @@
  * @param RecursionLevel contains current recursion depth
  * @param Chord baseline vector of first and second point
- * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
+ * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
  * @param OldNormal normal of the triangle which the baseline belongs to
+ * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere
  * @param Opt_Candidate candidate reference to return
- * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
  * @param Storage array containing two angles of current Opt_Candidate
  * @param RADIUS radius of ball
  * @param mol molecule structure with atoms and bonds
  */
+
+void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
+    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
+{
+
+  cout << Verbose(1) << "Candidate is "<< Candidate->nr << endl;
+  cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
+  /* OldNormal is normal vector on the old triangle
+   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck;
+  Vector TempNormal, Umkreismittelpunkt;
+  Vector Mittelpunkt;
+
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+  double BallAngle;
+  atom *Walker; // variable atom point
+
+
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  AngleCheck.CopyVector(&(Candidate->x));
+  AngleCheck.SubtractVector(&(a->x));
+  AngleCheck.ProjectOntoPlane(Chord);
+
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+
+  if (a != Candidate and b != Candidate)
+    {
+
+      Umkreisradius = SideA / 2.0 / sin(alpha);
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
+
+      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
+        {
+          sign = AngleCheck.ScalarProduct(direction1);
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+
+          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, sign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+          AngleCheck.CopyVector(&ReferencePoint);
+          AngleCheck.Scale(-1);
+          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          AngleCheck.AddVector(&Mittelpunkt);
+          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+          BallAngle = AngleCheck.Angle(OldNormal);
+          cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          sign = AngleCheck.ScalarProduct(direction1);
+          sign /= fabs(sign);
+
+
+          if (sign >0)
+            {
+              if (Storage[0]< -1.5) // first Candidate at all
+                {
+
+                  cout << "Next better candidate is " << *Candidate << " with ";
+                  Opt_Candidate = Candidate;
+                  Storage[0] = sign;
+                  Storage[1] = BallAngle;
+                  cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                    << Storage[0] << endl;
+
+
+                }
+              else
+                {
+                  if ( Storage[1] > BallAngle)
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = BallAngle;
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                        << Storage[0] << endl;
+                    }
+                  else
+                    {
+                      if (DEBUG)
+                        {
+                          cout << "Looses to better candidate" << endl;
+                        }
+                    }
+                }
+              }
+            else
+              {
+                cout << "Refused due to sign which is " << sign << endl;
+              }
+          }
+        else
+          {
+            if (DEBUG)
+              {
+                cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+              }
+          }
+      }
+    else
+      {
+        if (DEBUG)
+          {
+            cout << "identisch mit Ursprungslinie" << endl;
+          }
+      }
+
+
+
+  if (RecursionLevel < 7) // Five is the recursion level threshold.
+    {
+      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
+        {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
+              Candidate);
+          if (Walker == Parent)
+            { // don't go back the same bond
+              continue;
+            }
+          else
+            {
+              Find_next_suitable_point_via_Angle_of_Sphere(a, b, Walker, Candidate, RecursionLevel
+                  + 1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS,
+                  mol); //call function again
+            }
+        }
+    }
+}
+;
+
+
+  /** This recursive function finds a third point, to form a triangle with two given ones.
+   * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
+   *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
+   *  upon which we operate.
+   *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
+   *  direction and angle into Storage.
+   *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
+   *  with all neighbours of the candidate.
+   * @param a first point
+   * @param b second point
+   * @param Candidate base point along whose bonds to start looking from
+   * @param Parent point to avoid during search as its wrong direction
+   * @param RecursionLevel contains current recursion depth
+   * @param Chord baseline vector of first and second point
+   * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
+   * @param OldNormal normal of the triangle which the baseline belongs to
+   * @param Opt_Candidate candidate reference to return
+   * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
+   * @param Storage array containing two angles of current Opt_Candidate
+   * @param RADIUS radius of ball
+   * @param mol molecule structure with atoms and bonds
+   */
+
 void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
@@ -1669 +1884 @@
                   + 1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS,
                   mol); //call function again
+
             }
         }
@@ -1684 +1900 @@
  * @param N number of found triangles
  */
-void
-Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
+void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
     molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
     const double& RADIUS, int N)
@@ -1727 +1942 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-  {
-    Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
-    double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-
-    SideA = dif_b.Norm();
-      SideB = dif_a.Norm();
-      SideC = Chord->Norm();
-      //Chord->Scale(-1);
-
-      alpha = Chord->Angle(&dif_a);
-      beta = M_PI - Chord->Angle(&dif_b);
-      gamma = dif_a.Angle(&dif_b);
-
-      Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-      Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-  }
+
+    Vector Umkreismittelpunkt, a, b, c;
+    double alpha, beta, gamma;
+    a.CopyVector(&(T.endpoints[0]->node->x));
+    b.CopyVector(&(T.endpoints[1]->node->x));
+    c.CopyVector(&(T.endpoints[2]->node->x));
+    a.SubtractVector(&(T.endpoints[1]->node->x));
+    b.SubtractVector(&(T.endpoints[2]->node->x));
+    c.SubtractVector(&(T.endpoints[0]->node->x));
+
+    alpha = M_PI - a.Angle(&c);
+    beta = M_PI - b.Angle(&a);
+    gamma = M_PI - c.Angle(&b);
+
+    Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+
+
   cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
 
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
       Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
-      &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
+      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
       Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
-      &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
 
 
@@ -1799 +2018 @@
 ;
 
-void
-Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[2],
     molecule* mol, double RADIUS)
@@ -1860 +2078 @@
 ;
 
-void
-Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
+void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
   cout << Verbose(1) << "Looking for starting triangle \n";
@@ -1873 +2090 @@
   Vector helper;
   Vector Chord;
-  Vector Opt_Mittelpunkt;
+  Vector CenterOfFirstLine;
 
   Oben.Zero();
@@ -1936 +2153 @@
   // look in one direction of baseline for initial candidate
   Opt_Candidate = NULL;
-  Opt_Mittelpunkt;
-  Find_next_suitable_point(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
-      &Chord, &helper, &Oben, Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+  CenterOfFirstLine.CopyVector(&Chord);
+  CenterOfFirstLine.Scale(-0.5);
+  CenterOfFirstLine.AddVector(&(SecondPoint->x));
+
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
+      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
   // look in other direction of baseline for possible better candidate
-  Find_next_suitable_point(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
-      &Chord, &helper, &Oben, Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
+      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
   cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
 
@@ -1967 +2187 @@
 ;
 
-void
-Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
 {
   int N = 0;