Changeset 1ffa21


Ignore:
Timestamp:
Nov 27, 2008, 9:59:27 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:
69eb71
Parents:
03648b
Message:

The border.cpp contains the same functions as the changes in boundary, in case it is desired to have the convex and non convex routine seperated at a later time.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/border.cpp

    r03648b r1ffa21  
    11#include "molecules.hpp"
    22#include "boundary.hpp"
    3 
    4 
    5 
    6 
    73
    84
    … …  
    128   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
    139   */
    14   Vector *dif_a;
    15   Vector *dif_b;
    16   Vector *Chord;
    17   Vector *AngleCheck;
     10  Vector dif_a;
     11  Vector dif_b;
     12  Vector Chord;
     13  Vector AngleCheck;
    1814  atom *Walker;
    1915
    20   dif_a.CopyVector(a.x);
    21   dif_a.SubtractVector(Candidate->x);
    22   dif_b.CopyVector(b.x);
    23   dif.b.SubtractVector(Candidate->x);
    24   Chord.CopyVector(a.x);
    25   Chord.SubtractVector(b.x);
    26   AngleCheck.CopyVector(dif_a);
    27   AngleCheck.ProjectOntoPlane(Chord);
     16  dif_a.CopyVector(&a.x);
     17  dif_a.SubtractVector(&Candidate.x);
     18  dif_b.CopyVector(&b.x);
     19  dif_b.SubtractVector(&Candidate.x);
     20  Chord.CopyVector(&a.x);
     21  Chord.SubtractVector(&b.x);
     22  AngleCheck.CopyVector(&dif_a);
     23  AngleCheck.ProjectOntoPlane(&Chord);
    2824
    2925  //Storage eintrag fuer aktuelles Atom
    30   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
     26  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
    3127  {
    3228
    … …  
    4440          {
    4541            Storage[0]=(double)Candidate.nr;
    46             Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif.ScalarProduct(d1));
     42            Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
    4743            Storage[2]=AngleCheck.Angle(d2);
    4844          }
    … …  
    5248  if (n<5)
    5349    {
    54       for(i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
    55         {
    56           while (Candidate.nr != mol->ListOfBonds[Candidate.nr][i])
     50      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
     51        {
     52          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))
    5753            {
    58               Walker = Walker.next;
     54              Walker = Walker->next;
    5955            }
    6056
    61           Find_next_suitable_point(a, b, Walker, n+1, d1, d2, RADIUS, mol);
    62         }
    63     }
    64 }
    65 
    66 
    67 void Find_next_suitable_triangle(molecule mol, BoundaryLineSet Line, BoundaryTriangleSet T)
    68 {
    69   Vector CenterOfLine = Line->endpoints.node[0].x;
     57          Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
     58        }
     59    }
     60};
     61
     62
     63void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
     64{
     65  Vector CenterOfLine = Line.endpoints[0]->node->x;
    7066  Vector direction1;
    7167  Vector direction2;
    7268  Vector helper;
    73 
    74   double *Storage[3];
    75   Storage[0]=-1;   // Id must be positive, we see should nothing be done
    76   Storage[1]=-1;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
    77   Storage[2]=-10;  // This is also lower then any value produced by an eligible atom, which are all positive
    78 
    79  
    80   helper.CopyVector(Line->endpoints[0].x);
     69  atom* Walker;
     70
     71  double Storage[3];
     72  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
     73  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
     74  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
     75
     76 
     77  helper.CopyVector(&(Line.endpoints[0]->node->x));
    8178  for (int i =0; i<3; i++)
    8279    {
    83       if (T->endpoints[i].node.nr != Line->endpoints[0].node.nr && T->endpoints[i].node.nr!=Line->endpoints[1].node.nr)
    84         {
    85           helper.SubtractVector(T->endpoints[i].x);
     80      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
     81        {
     82          helper.SubtractVector(&T.endpoints[i]->node->x);
    8683          break;
    8784        }
    … …  
    8986
    9087
    91   direction1.CopyVector(Line->endpoints.node[0].x);
    92   direction1.Subtract(Line->endpoints.node[1].x);
    93   direction1.Crossproduct(T.NormalVector);
    94 
    95   if (direction1.ScalarProduct(helper)>0)
     88  direction1.CopyVector(&Line.endpoints[0]->node->x);
     89  direction1.SubtractVector(&Line.endpoints[1]->node->x);
     90  direction1.VectorProduct(T.NormalVector);
     91
     92  if (direction1.ScalarProduct(&helper)>0)
    9693    {
    9794      direction1.Scale(-1);
    9895    }
    9996
    100   Find_next_suitable_point(Line->endpoints.node[0], Line->endpoints.node[1], Line->endpoints->node[0], 0, direction1, T.NormalVector, Storage, RADIUS);
     97  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
    10198 
    10299  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
    103100  // Next Triangle is Line, atom with number in Storage[0]
    104101
    105   Walker=mol->start;
    106   while (Walker.nr != (int)Storage[0])
     102  Walker= mol->start;
     103  while (Walker->nr != (int)Storage[0])
    107104    {
    108105      Walker = Walker->next;
    … …  
    111108  AddPoint(Walker);
    112109
    113   BPS[0] = BoundaryPointSet(Walker);
    114   BPS[1] = Line->enpoints.node[0];
     110  BPS[0] = new class BoundaryPointSet(Walker);
     111  BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
    115112  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    116   BPS[0] = Walker;
    117   BPS[1] = Line->endpoints.node[1];
     113  BPS[0] = new class BoundaryPointSet(Walker);
     114  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
    118115  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    119   BLS[2] = line;
     116  BLS[2] = new class BoundaryLineSet(Line);
    120117
    121118  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    … …  
    125122  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
    126123    {
    127       if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary->end)
     124      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
    128125        {
    129126          LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
    … …  
    131128        }
    132129    }
    133   GetNormalVector(BTS.NormalVector);
    134 
    135   if( (BTS.NormalVector.ScalarProduct(T.NormalVecotr)<0 && Storage[1]>0) || \
    136       (BTS.NormalVector.ScalarProduct(T.NormalVecotr)>0 && Storage[1]<0))
    137     {
    138       BTS.NormalVector.Scale(-1);
    139     }
    140  
    141 }
    142 
    143 
    144 void Find_starting_triangle(molecule mol)
    145 {
     130  BTS->GetNormalVector(*BTS->NormalVector);
     131
     132  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
     133      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
     134    {
     135      BTS->NormalVector->Scale(-1);
     136    }
     137 
     138};
     139
     140
     141void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
     142{
     143  int i;
     144  Vector *AngleCheck;
     145  atom* Walker;
     146
     147  AngleCheck->CopyVector(&Candidate.x);
     148  AngleCheck->SubtractVector(&a.x);
     149  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
     150    {
     151      Storage[0]=(double)(Candidate.nr);
     152      Storage[1]=AngleCheck->ScalarProduct(&Oben);
     153    };
     154
     155  if (n<5)
     156    {
     157      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
     158        {
     159          Walker = mol.start;
     160          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))
     161            {
     162              Walker = Walker->next;
     163            };
     164         
     165          Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
     166            };
     167    };
     168 
     169
     170};
     171
     172
     173void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
     174{
     175  int i=0;
    146176  atom Walker;
    147177  atom Walker2;
    … …  
    152182  Vector helper;
    153183
    154   Oben.Zero;
    155 
    156 
    157   for(int i =0; i<3; i++)
     184  Oben.Zero();
     185
     186
     187  for(i =0; i<3; i++)
    158188    {
    159189      max_index[i] =-1;
    … …  
    161191    }
    162192
    163   Walker = mol->start;
    164   while (Walker->next != NULL)
     193  Walker = *mol.start;
     194  while (Walker.next != NULL)
    165195    {
    166196      for (i=0; i<3; i++)
    167197        {
    168           if (Walker.x[i]>max_coordinate[i])
     198          if (Walker.x.x[i]>max_coordinate[i])
    169199            {
    170               max_coordinate[i]=Walker.x[i];
     200              max_coordinate[i]=Walker.x.x[i];
    171201              max_index[i]=Walker.nr;
    172202            }
    … …  
    178208
    179209  Oben.x[k]=1;
    180   Walker = mol->start;
     210  Walker = *mol.start;
    181211  while (Walker.nr != max_index[k])
    182212    {
    183       Walker =Walker->next;
    184     }
    185 
    186   double *Storage[3];
    187   Storage[0]=-1;   // Id must be positive, we see should nothing be done
    188   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.
    189   Storage[2]=-10;  // This will be an angle looking for the third point.
    190 
    191 
    192   for (i=0; i< mol->NumberOfBondsPerAtoms[Walker.nr]; i++)
    193     {
    194       Walker2 = mol->start;
    195       while (Walker2.nr != mol->ListOfBondsPerAtoms[Walker.nr][i])  // Stimmt die Ueberpruefung $$$
    196         {
    197           Walker2 =Walker2->next;
    198         }
    199 
    200       Find_second_point_for_Tesselation(Walker, Walker2, Oben, Storage);
    201     }
    202 
    203   Walker2=mol->start;
     213      Walker = *Walker.next;
     214    }
     215
     216  double Storage[3];
     217  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
     218  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.
     219  Storage[2]=-10.;  // This will be an angle looking for the third point.
     220
     221
     222  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
     223    {
     224      Walker2 = *mol.start;
     225      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) )
     226        {
     227          Walker2 = *Walker2.next;
     228        }
     229
     230      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
     231    }
     232
     233  Walker2 = *mol.start;
    204234
    205235  while (Walker2.nr != int(Storage[0]))
    206236    {
    207       Walker = Walker.next;
    208     }
    209  
    210   helper.copyVector(Walker.x);
    211   helper.Subtract(Walker2.x);
    212   Oben.ProjectOntoPlane(helper.x);
    213   helper.VectorProduct(Oben);
    214 
    215   Find_next_suitable_point(Walker, Walker2, Candidate, 0, helper, Oben, Storage, Radius);
    216   Walker3 = mol->start;
     237      Walker = *Walker.next;
     238    }
     239 
     240  helper.CopyVector(&Walker.x);
     241  helper.SubtractVector(&Walker2.x);
     242  Oben.ProjectOntoPlane(&helper);
     243  helper.VectorProduct(&Oben);
     244
     245       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);
     246  Walker3 = *mol.start;
    217247  while (Walker3.nr != int(Storage[0]))
    218248    {
    219       Walker3 = Walker3->next;
     249      Walker3 = *Walker3.next;
    220250    }
    221251
    222252  //Starting Triangle is Walker, Walker2, index Storage[0]
    223253
    224   AddPoint(Walker);
    225   AddPoint(Walker2);
    226   AddPoint(Walker3);
    227 
    228   BPS[0] = BoundaryPointSet(Walker);
    229   BPS[1] = BoundaryPointSet(Walker2);
     254  AddPoint(&Walker);
     255  AddPoint(&Walker2);
     256  AddPoint(&Walker3);
     257
     258  BPS[0] = new class BoundaryPointSet(&Walker);
     259  BPS[1] = new class BoundaryPointSet(&Walker2);
    230260  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    231   BPS[0] = Walker;
    232   BPS[1] = Walker3;
     261  BPS[0] = new class BoundaryPointSet(&Walker);
     262  BPS[1] = new class BoundaryPointSet(&Walker3);
    233263  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    234   BPS[0] = Walker;
    235   BPS[1] = Walker2;
     264  BPS[0] = new class BoundaryPointSet(&Walker);
     265  BPS[1] = new class BoundaryPointSet(&Walker2);
    236266  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount); 
    237267
    … …  
    244274      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
    245275      LinesOnBoundaryCount++;
    246     }
    247 
    248 }
    249 
    250 
    251 void Find_non_convex_border(Tesselation Tess, molecule mol) // this needs input of type molecule
    252 {
    253 
    254   Tess.Find_starting_triangle(mol);
    255 
    256   for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
     276    };
     277
     278       BTS->GetNormalVector(*BTS->NormalVector);
     279
     280       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
     281         {
     282           BTS->NormalVector->Scale(-1);
     283         }
     284};
     285
     286
     287void Find_non_convex_border(Tesselation* Tess, molecule mol)
     288{
     289  const double RADIUS =6;
     290  Tess->Find_starting_triangle(mol, RADIUS);
     291
     292  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++)
    257293    if (baseline->second->TrianglesCount == 1)
    258294      {
    259         Find_next_suitable_triangle(mol, baseline->second, baseline->second->triangles.begin()->second); //the line is there, so there is a triangle, but only one.
     295        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.
    260296
    261297      }
    … …  
    264300        printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
    265301      }
    266 }
     302};
Note: See TracChangeset for help on using the changeset viewer.