Changeset 03648b for src/boundary.cpp

Timestamp:

Nov 27, 2008, 9:55:08 AM (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:

1ffa21

Parents:

51695c

Message:

In vector a function for calculation of the vector-(cross-)product has been added.
In Boundary a new way for finding the non-convex boundary is implemented.
Currently problem with comparison of the return value of the map::find routine.

File:

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

src/boundary.cpp

@@ -1 +1 @@
 #include "molecules.hpp"
 #include "boundary.hpp"
+
+
+// =================spaeter gebrauchte
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol);
+void Find_non_convex_border(Tesselation Tess, molecule mol);
+
 
 // ======================================== Points on Boundary =================================
@@ -471 +477 @@
   double volume = 0.;
   double PyramidVolume = 0.;
-  double G,h;
+  double G,h; 
   Vector x,y;
   double a,b,c;
 
+  Find_non_convex_border(*TesselStruct, *mol);
+  /*
   // 1. calculate center of gravity
   *out << endl;
@@ -518 +526 @@
 
   *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
+  */
+
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
@@ -539 +549 @@
   *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 
-
+  /*
   // 7. translate all points back from CoG
   *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
@@ -548 +558 @@
     Walker->x.Translate(CenterOfGravity);
   }
-  
+  */
+
+
+
+
+
   // 8. Store triangles in tecplot file
   if (tecplot != NULL) {
@@ -994 +1009 @@
     PointsOnBoundaryCount++;
 };
+
+
+
+
+
+
+//====================================================================================================================
+//
+// ab hier habe ich das verzapft.
+//
+//====================================================================================================================
+
+
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol)
+{
+  /* d2 ist der Normalenvektor auf dem Dreieck,
+   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
+   */
+  Vector dif_a;
+  Vector dif_b;
+  Vector Chord;
+  Vector AngleCheck;
+  atom *Walker;
+
+  dif_a.CopyVector(&a.x);
+  dif_a.SubtractVector(&Candidate.x);
+  dif_b.CopyVector(&b.x);
+  dif_b.SubtractVector(&Candidate.x);
+  Chord.CopyVector(&a.x);
+  Chord.SubtractVector(&b.x);
+  AngleCheck.CopyVector(&dif_a);
+  AngleCheck.ProjectOntoPlane(&Chord);
+
+  //Storage eintrag fuer aktuelles Atom
+  if (Chord.Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  {
+
+    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
+      {
+        Storage[0]=(double)Candidate.nr;
+        Storage[1]=1;
+        Storage[2]=AngleCheck.Angle(d2);
+      }
+    else
+      {
+        if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
+            (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
+          //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+          {
+            Storage[0]=(double)Candidate.nr;
+            Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+            Storage[2]=AngleCheck.Angle(d2);
+          }
+      }
+  }
+
+  if (n<5)
+    {
+      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+        {
+          while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+            {
+              Walker = Walker->next;
+            }
+
+          Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
+        }
+    }
+};
+
+
+void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
+{
+  Vector CenterOfLine = Line.endpoints[0]->node->x;
+  Vector direction1;
+  Vector direction2;
+  Vector helper;
+  atom* Walker;
+
+  double Storage[3];
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
+
+  
+  helper.CopyVector(&(Line.endpoints[0]->node->x));
+  for (int i =0; i<3; i++)
+    {
+      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
+        {
+          helper.SubtractVector(&T.endpoints[i]->node->x);
+          break;
+        }
+    }
+
+
+  direction1.CopyVector(&Line.endpoints[0]->node->x);
+  direction1.SubtractVector(&Line.endpoints[1]->node->x);
+  direction1.VectorProduct(T.NormalVector);
+
+  if (direction1.ScalarProduct(&helper)>0)
+    {
+      direction1.Scale(-1);
+    }
+
+  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
+  
+  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
+  // Next Triangle is Line, atom with number in Storage[0]
+
+  Walker= mol->start;
+  while (Walker->nr != (int)Storage[0])
+    {
+      Walker = Walker->next;
+    }
+
+  AddPoint(Walker);
+
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BLS[2] = new class BoundaryLineSet(Line);
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
+    {
+      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
+        {
+          LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+          LinesOnBoundaryCount++;
+        }
+    }
+  BTS->GetNormalVector(*BTS->NormalVector);
+
+  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
+      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+    {
+      BTS->NormalVector->Scale(-1);
+    }
+  
+};
+
+
+void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
+{
+  int i;
+  Vector *AngleCheck;
+  atom* Walker;
+
+  AngleCheck->CopyVector(&Candidate.x);
+  AngleCheck->SubtractVector(&a.x);
+  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
+    {
+      Storage[0]=(double)(Candidate.nr);
+      Storage[1]=AngleCheck->ScalarProduct(&Oben);
+    };
+
+  if (n<5)
+    {
+      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
+        {
+          Walker = mol.start;
+          while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+            {
+              Walker = Walker->next;
+            };
+          
+          Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
+            };
+    };
+  
+
+};
+
+
+void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
+{
+  int i=0;
+  atom Walker;
+  atom Walker2;
+  atom Walker3;
+  int max_index[3];
+  double max_coordinate[3];
+  Vector Oben;
+  Vector helper;
+
+  Oben.Zero();
+
+
+  for(i =0; i<3; i++)
+    {
+      max_index[i] =-1;
+      max_coordinate[i] =-1;
+    }
+
+  Walker = *mol.start;
+  while (Walker.next != NULL)
+    {
+      for (i=0; i<3; i++)
+        {
+          if (Walker.x.x[i]>max_coordinate[i])
+            {
+              max_coordinate[i]=Walker.x.x[i];
+              max_index[i]=Walker.nr;
+            }
+        }
+    }
+
+  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
+  const int k=0;
+
+  Oben.x[k]=1;
+  Walker = *mol.start;
+  while (Walker.nr != max_index[k])
+    {
+      Walker = *Walker.next;
+    }
+
+  double Storage[3];
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[2]=-10.;  // This will be an angle looking for the third point.
+
+
+  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
+    {
+      Walker2 = *mol.start;
+      while (Walker2.nr != (mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom->nr : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr) ) 
+        {
+          Walker2 = *Walker2.next;
+        }
+
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
+    }
+
+  Walker2 = *mol.start;
+
+  while (Walker2.nr != int(Storage[0]))
+    {
+      Walker = *Walker.next;
+    }
+  
+  helper.CopyVector(&Walker.x);
+  helper.SubtractVector(&Walker2.x);
+  Oben.ProjectOntoPlane(&helper);
+  helper.VectorProduct(&Oben);
+
+       Find_next_suitable_point(Walker, Walker2, *(mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom), 0, &helper, &Oben, Storage, RADIUS, mol);
+  Walker3 = *mol.start;
+  while (Walker3.nr != int(Storage[0]))
+    {
+      Walker3 = *Walker3.next;
+    }
+
+  //Starting Triangle is Walker, Walker2, index Storage[0]
+
+  AddPoint(&Walker);
+  AddPoint(&Walker2);
+  AddPoint(&Walker3);
+
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker2);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker3);
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker2);
+  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) 
+    {
+      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+      LinesOnBoundaryCount++;
+    };
+
+       BTS->GetNormalVector(*BTS->NormalVector);
+
+       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+         {
+           BTS->NormalVector->Scale(-1);
+         }
+};
+
+
+void Find_non_convex_border(Tesselation* Tess, molecule mol) 
+{
+  const double RADIUS =6;
+  Tess->Find_starting_triangle(mol, RADIUS);
+
+  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++) 
+    if (baseline->second->TrianglesCount == 1) 
+      {
+        Tess->Find_next_suitable_triangle(&mol, baseline->second, baseline->second->triangles.begin()->second, RADIUS); //the line is there, so there is a triangle, but only one.
+
+      }
+    else
+      {
+        printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+      }
+};