Context Navigation

-              ra8bcea6
+              rf714979
   double a,b,c;
   Find_non_convex_border(out, tecplot, *TesselStruct, mol);
+  /*
+  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
   // 1. calculate center of gravity
   *out << endl;
 …
   *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
+  */
 …
   *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;
 …
     Walker->x.Translate(CenterOfGravity);
+  }
+  */
 …
  *  with all neighbours of the candidate.
  */
 void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
+{
   /* d2 is normal vector on the triangle
    * d1 is normal on the triangle line, from which we come, as well as on d2.
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *OldNormal, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
+{
+  /* OldNormal is normal vector on the old triangle
+   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
    */
   Vector dif_a; //Vector from a to candidate
 …
+  {
           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 (Chord->Norm()/(2*sin(0.5*dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find if Ball will touch atom
+          {
                   if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
 …
                           Opt_Candidate = Candidate;
                           Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
                           Storage[1]=AngleCheck.Angle(d2);
+                          Storage[1]=AngleCheck.Angle(OldNormal);
+                  }
                   else
+                  {
                           if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]< AngleCheck.Angle(d2)) or \
                                           (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]> AngleCheck.Angle(d2)))
+                          if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]> AngleCheck.Angle(OldNormal)) or \
+                                          (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]< AngleCheck.Angle(OldNormal)))
                                   //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+                          {
                                   Opt_Candidate = Candidate;
                                   Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
                                   Storage[1]=AngleCheck.Angle(d2);
+                                  Storage[1]=AngleCheck.Angle(OldNormal);
+                          }
                           else
+                          {
                                   if (problem)
                                           cout << "Loses to better candidate" << endl;
+                                          cout << "Looses to better candidate" << endl;
+                          }
+                  }
 …
                   Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
                   Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
+                  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, OldNormal, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
+      }
+  }
 …
  * this function fins a triangle to a line, adjacent to an existing one.
  */
 void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS)
+void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N)
+{
         cout << Verbose(1) << "Looking for next suitable triangle \n";
 …
   direction1.VectorProduct(&(T.NormalVector));
   if (direction1.ScalarProduct(&helper)>0)
+  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, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
   if (Opt_Candidate==NULL)
+  if (N>-1)
+  {
           cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
 …
   cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 …
   cout << Verbose(1) << "Checking for already present lines ... " << endl;
+  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
+    {
+  if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
+        {
+                LinesOnBoundaryCount++;
+        }
+    }
+  for(int i=0; i<NDIM; i++) // sind Linien bereits vorhanden ???
+  {
+          if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
+          {
+                  LinesOnBoundaryCount++;
+          }
+  }
   cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
 …
         atom* Walker;
+        AngleCheck.CopyVector(&(Candidate->x));
+        AngleCheck.SubtractVector(&(a->x));
+        if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
+    {
+          //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Opt_Candidate=Candidate;
+      Storage[0]=AngleCheck.Angle(&Oben);
+      //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+    };
+        if (a->nr !=Candidate->nr)
+        {
+                AngleCheck.CopyVector(&(Candidate->x));
+                AngleCheck.SubtractVector(&(a->x));
+                if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
+                {
+                        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                        Opt_Candidate=Candidate;
+                        Storage[0]=AngleCheck.Angle(&Oben);
+                        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+                }
+                else{
+                        if (AngleCheck.Norm() > RADIUS)
+                        {
+                                cout << Verbose(1) << "refused due to Radius" << AngleCheck.Norm() << endl;
+                        }
+                        else{
+                                cout << Verbose(1) << "Supposedly looses to better candidate" << Opt_Candidate->nr << endl;
+                        }
+                }
+        }
         if (n<5)
+    {
 …
   int i=0;
   atom* Walker;
   atom* Walker2;
   atom* Walker3;
+  atom* FirstPoint;
+  atom* SecondPoint;
   int max_index[3];
   double max_coordinate[3];
 …
+    {
         Walker = Walker->next;
       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;
+            }
+      }
+    }
 cout << Verbose(1) << "Found starting atom \n";
+        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;
+                }
+        }
+    }
+  cout << Verbose(1) << "Found maximum coordinates. "<< endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
   const int k=2;
+  const int k=0;
   Oben.x[k]=1.;
   Walker = mol->start;
   Walker = Walker->next;
   while (Walker->nr != max_index[k])
+  FirstPoint = mol->start;
+  FirstPoint = FirstPoint->next;
+  while (FirstPoint->nr != max_index[k])
+    {
       Walker = Walker->next;
+    }
 cout << Verbose(1) << "Number of start atom: " << Walker->nr << endl;
+          FirstPoint = FirstPoint->next;
+    }
+  cout << Verbose(1) << "Coordinates of start atom: " << FirstPoint->x.x[0] << endl;
   double Storage[2];
   atom* Opt_Candidate = NULL;
   Storage[0]=999999.;   // 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[1]=999999.;  // This will be an angle looking for the third point.
+cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[Walker->nr] << endl;
+  for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
+    {
+          Walker2 = mol->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+  Walker2 = Opt_Candidate;
+  cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
+  Find_second_point_for_Tesselation(FirstPoint, mol->ListOfBondsPerAtom[FirstPoint->nr][0]->GetOtherAtom(FirstPoint), 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
+  SecondPoint = Opt_Candidate;
   Opt_Candidate=NULL;
   helper.CopyVector(&(Walker->x));
   helper.SubtractVector(&(Walker2->x));
+  helper.CopyVector(&(FirstPoint->x));
+  helper.SubtractVector(&(SecondPoint->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
 …
   Storage[1]= 9999999.; // This will be an angle looking for the third point.
   Chord.CopyVector(&(Walker->x)); // bring into calling function
   Chord.SubtractVector(&(Walker2->x));
+  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+  Chord.SubtractVector(&(SecondPoint->x));
   cout << Verbose(1) << "Looking for third point candidates \n";
   Find_next_suitable_point(Walker, Walker2, mol->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
   //Starting Triangle is Walker, Walker2, Opt_Candidate
   AddPoint(Walker);
   AddPoint(Walker2);
+  Find_next_suitable_point(FirstPoint, SecondPoint, mol->ListOfBondsPerAtom[FirstPoint->nr][0]->GetOtherAtom(FirstPoint), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
+  //Starting Triangle is Walker, SecondPoint, Opt_Candidate
+  AddPoint(FirstPoint);
+  AddPoint(SecondPoint);
   AddPoint(Opt_Candidate);
   cout << Verbose(1) << "The found starting triangle consists of " << *Walker << ", " << *Walker2 << " and " << *Opt_Candidate << "." << endl;
   BPS[0] = new class BoundaryPointSet(Walker);
   BPS[1] = new class BoundaryPointSet(Walker2);
+  cout << Verbose(1) << "The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << "." << endl;
+  BPS[0] = new class BoundaryPointSet(FirstPoint);
+  BPS[1] = new class BoundaryPointSet(SecondPoint);
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[0] = new class BoundaryPointSet(FirstPoint);
   BPS[1] = new class BoundaryPointSet(Opt_Candidate);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
   BPS[1] = new class BoundaryPointSet(Walker2);
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  BPS[1] = new class BoundaryPointSet(SecondPoint);
   BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 …
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol)
+{
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
+{
+        int N =0;
+        struct Tesselation *Tess = new Tesselation;
         cout << Verbose(1) << "Entering finding of non convex hull. " << endl;
         cout << flush;
+  const double RADIUS =6;
+  LineMap::iterator baseline;
+  Tess.Find_starting_triangle(mol, RADIUS);
+  baseline = Tess.LinesOnBoundary.begin();
+  while (baseline != Tess.LinesOnBoundary.end()) {
+    if (baseline->second->TrianglesCount == 1)
+      {
+        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+                Tess.Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+        cout << Verbose(1) << "End of Tesselation ... " << endl;
+      }
+    else
+      {
+        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
+      }
+  baseline++;
+  }
+};
+        const double RADIUS =6.;
+        LineMap::iterator baseline;
+        Tess->Find_starting_triangle(mol, RADIUS);
+        baseline = Tess->LinesOnBoundary.begin();
+        while (baseline != Tess->LinesOnBoundary.end()) {
+                if (baseline->second->TrianglesCount == 1)
+                {
+                        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+                        Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
+                        cout << Verbose(1) << "End of Tesselation ... " << endl;
+                }
+                else
+                {
+                        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
+                }
+                N++;
+                baseline++;
+        }
+};

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset f714979 for src/boundary.cpp

Legend:

src/boundary.cpp

Download in other formats: