Changeset 10af0d

Timestamp:

Dec 17, 2008, 3:56:07 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:

e9fa06f

Parents:

e4ea46

Message:

Solved possible problem with Second quadrant overhaul

File:

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

src/boundary.cpp

@@ -1444 +1444 @@
   /* OldNormal is normal vector on the old triangle
    * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
-   * Chord points from a to b.
+   * Chord points from b to a!!!
    */
   Vector dif_a; //Vector from a to candidate
@@ -1495 +1495 @@
 
       Umkreisradius = SideA / 2.0 / sin(alpha);
-      cout << Umkreisradius << endl;
-      cout << SideB / 2.0 / sin(beta) << endl;
-      cout << SideC / 2.0 / sin(gamma) << endl;
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
 
       if (Umkreisradius < RADIUS)
@@ -1510 +1510 @@
           Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
 
-          cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-          cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-          cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-          cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
 
           TempNormal.CopyVector(&dif_a);
@@ -1523 +1519 @@
           TempNormal.Normalize();
           Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-          cout << "Umkreisradius ist " << Umkreisradius << endl;
-          cout << "Restradius ist " << Restradius << endl;
+
           TempNormal.Scale(Restradius);
-          cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+
           Mittelpunkt.CopyVector(&Umkreismittelpunkt);
           Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-          cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-          cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-          cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-          cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+
+
 
           if (Storage[0]< -1.5) // first Candidate at all
@@ -1547 +1537 @@
               cout << "Angle is " << Storage[1] << ", Halbraum ist "
                   << Storage[0] << endl;
+
+
             }
           else
@@ -1554 +1546 @@
 
               Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-              cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+              //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
 
 
               if (Distance >RADIUS) // We have no interference and may now check whether the new point is better.
                 {
-                  cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
-
-                  if ((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))// || //This will give absolute preference to those in "right-hand" quadrants
-                      //sqrt(Candidate->x.Distance(Opt_Mittelpunkt)) < RADIUS)    //and those where Candidate would be within old Sphere.
+                  //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+
+                  if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon)) || //This will give absolute preference to those in "right-hand" quadrants
+                      (Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
                     {
                       cout << "Next better candidate is " << *Candidate << " with ";
@@ -1571 +1563 @@
                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
                           << Storage[0] << endl;
+/*
+                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
+                                    cout << "Restradius ist " << Restradius << endl;
+                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+                                    */
+
                     }
                   else
@@ -1592 +1602 @@
               else
                 {
-                  if (DEBUG)
-                    cout << "Looses to better candidate" << endl;
+                  if (sign>0 && AngleCheck.Angle(OldNormal)>0 && Storage[0]<0)
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                      << Storage[0] << endl;
+
+/*
+                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
+                                    cout << "Restradius ist " << Restradius << endl;
+                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+                                    */
+
+
+                    }
+                  else
+                    {
+                      if (DEBUG)
+                        cout << "Looses to better candidate" << endl;
+                    }
                 }
             }
@@ -1612 +1655 @@
     }
 
-  if (RecursionLevel < 7) // Five is the recursion level threshold.
+  if (RecursionLevel < 9) // Five is the recursion level threshold.
     {
       for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
@@ -1685 +1728 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-  cout << Verbose(1) << "Looking for third point candidates for triangle ... "
-      << endl;
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+
   Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
       Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
@@ -1693 +1736 @@
       Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+
 
   if ((TrianglesOnBoundaryCount % 10) == 0)
@@ -1704 +1748 @@
   }
 
-  if ((Opt_Candidate == NULL) || (N<12))
+  if (TrianglesOnBoundaryCount >1000 )
     {
       cout << Verbose(1)
           << "No new Atom found, triangle construction will crash" << endl;
       write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
-      Find_next_suitable_point(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, Line.endpoints[1]->node,
-          Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector),
-          Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+      cout << "This is currently added candidate" << Opt_Candidate << endl;
     }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
@@ -1733 +1770 @@
   AddTriangleToLines();
   cout << "New triangle with " << *BTS << endl;
-
-  cout << Verbose(1) << "Constructing normal vector for this triangle ... "
-      << endl;
+  cout << "We have "<< TrianglesOnBoundaryCount << endl;
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
 
   BTS->GetNormalVector(BTS->NormalVector);
 
-  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)
-      || ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0]
-          < 0)))
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
+      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
     {
       BTS->NormalVector.Scale(-1);