Changeset 0d4daf for src/Fragmentation


Ignore:
Timestamp:
Aug 20, 2014, 1:04:08 PM (11 years ago)
Author:
Frederik Heber <heber@…>
Children:
90426a
Parents:
6ef48c
git-author:
Frederik Heber <heber@…> (05/30/14 13:31:15)
git-committer:
Frederik Heber <heber@…> (08/20/14 13:04:08)
Message:

Moved all code from header of SphericalPointDistribution into cpp module.

Location:
src/Fragmentation
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r6ef48c r0d4daf  
    1616
    1717#include "CodePatterns/Assert.hpp"
    18 #include "CodePatterns/IteratorAdaptors.hpp"
    19 #include "CodePatterns/toString.hpp"
    2018
    21 #include <algorithm>
    2219#include <cmath>
    23 #include <limits>
    2420#include <list>
    25 #include <vector>
    26 #include <map>
    2721
    28 #include "LinearAlgebra/Line.hpp"
    29 #include "LinearAlgebra/RealSpaceMatrix.hpp"
    3022#include "LinearAlgebra/Vector.hpp"
    3123
    … …  
    6557  }
    6658
     59
    6760  /** Matches a given spherical distribution with another containing more
    6861   * points.
    … …  
    9588};
    9689
    97 typedef std::list<unsigned int> IndexList_t;
    98 typedef std::vector<unsigned int> IndexArray_t;
    99 typedef std::vector<Vector> VectorArray_t;
    100 typedef std::vector<double> DistanceArray_t;
     90// declare specializations
    10191
    102 DistanceArray_t calculatePairwiseDistances(
    103     const std::vector<Vector> &_points,
    104     const IndexList_t &_indices
    105     )
    106 {
    107   DistanceArray_t result;
    108   for (IndexList_t::const_iterator firstiter = _indices.begin();
    109       firstiter != _indices.end(); ++firstiter) {
    110     for (IndexList_t::const_iterator seconditer = firstiter;
    111         seconditer != _indices.end(); ++seconditer) {
    112       if (firstiter == seconditer)
    113         continue;
    114       const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
    115       result.push_back(distance);
    116     }
    117   }
    118   return result;
    119 }
    120 
    121 // class generator: taken from www.cplusplus.com example std::generate
    122 struct c_unique {
    123   int current;
    124   c_unique() {current=0;}
    125   int operator()() {return ++current;}
    126 } UniqueNumber;
    127 
    128 /** Returns squared L2 error of the given \a _Matching.
    129  *
    130  * We compare the pair-wise distances of each associated matching
    131  * and check whether these distances each match between \a _old and
    132  * \a _new.
    133  *
    134  * \param _old first set of points (fewer or equal to \a _new)
    135  * \param _new second set of points
    136  * \param _Matching matching between the two sets
    137  * \return pair with L1 and squared L2 error
    138  */
    139 std::pair<double, double> calculateErrorOfMatching(
    140     const std::vector<Vector> &_old,
    141     const std::vector<Vector> &_new,
    142     const IndexList_t &_Matching)
    143 {
    144   std::pair<double, double> errors( std::make_pair( 0., 0. ) );
    145 
    146   if (_Matching.size() > 1) {
    147     // convert matching into two vectors to calculate distance among another
    148 
    149     // calculate all pair-wise distances
    150     IndexList_t keys(_Matching.size());
    151     std::generate (keys.begin(), keys.end(), UniqueNumber);
    152     const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
    153     const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
    154 
    155     ASSERT( firstdistances.size() == keys.size()*(keys.size()+1)/2,
    156         "calculateL2ErrorOfMatching() - unexpected pair-wise distance array size.");
    157     ASSERT( firstdistances.size() == seconddistances.size(),
    158         "calculateL2ErrorOfMatching() - mismatch in pair-wise distance array sizes.");
    159     DistanceArray_t::const_iterator firstiter = firstdistances.begin();
    160     DistanceArray_t::const_iterator seconditer = seconddistances.begin();
    161     for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
    162         ++firstiter, ++seconditer) {
    163       const double gap = *firstiter - *seconditer;
    164       // L1 error
    165       if (errors.first < gap)
    166         errors.first = gap;
    167       // L2 error
    168       errors.second += gap*gap;
    169     }
    170   }
    171 
    172   return errors;
    173 }
    174 
    175 SphericalPointDistribution::Polygon_t removeMatchingPoints(
    176     const SphericalPointDistribution::Polygon_t &_points,
    177     const IndexList_t &_matchingindices
    178     )
    179 {
    180   SphericalPointDistribution::Polygon_t remainingpoints;
    181   IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
    182   std::sort(indices.begin(), indices.end());
    183   IndexArray_t::const_iterator valueiter = indices.begin();
    184   SphericalPointDistribution::Polygon_t::const_iterator pointiter =
    185       _points.begin();
    186   for (unsigned int i=0; i< _points.size(); ++i, ++pointiter) {
    187     // skip all those in values
    188     if (*valueiter == i)
    189       ++valueiter;
    190     else
    191       remainingpoints.push_back(*pointiter);
    192   }
    193 
    194   return remainingpoints;
    195 }
    196 
    197 /** Rotates a given polygon around x, y, and z axis by the given angles.
    198  *
    199  * Essentially, we concentrate on the three points of the polygon to rotate
    200  * to the correct position. First, we rotate its center via \a angles,
    201  * then we rotate the "triangle" around itself/\a _RotationAxis by
    202  * \a _RotationAngle.
    203  *
    204  * \param _polygon polygon whose points to rotate
    205  * \param _angles vector with rotation angles for x,y,z axis
    206  * \param _RotationAxis
    207  * \param _RotationAngle
    208  */
    209 SphericalPointDistribution::Polygon_t rotatePolygon(
    210     const SphericalPointDistribution::Polygon_t &_polygon,
    211     const std::vector<double> &_angles,
    212     const Line &_RotationAxis,
    213     const double _RotationAngle)
    214 {
    215   SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
    216   RealSpaceMatrix rotation;
    217   ASSERT( _angles.size() == 3,
    218       "rotatePolygon() - not exactly "+toString(3)+" angles given.");
    219   rotation.setRotation(_angles[0] * M_PI/180., _angles[1] * M_PI/180., _angles[2] * M_PI/180.);
    220 
    221   // apply rotation angles
    222   for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
    223       iter != rotated_polygon.end(); ++iter) {
    224     *iter = rotation * (*iter);
    225     _RotationAxis.rotateVector(*iter, _RotationAngle);
    226   }
    227 
    228   return rotated_polygon;
    229 }
    230 
    231 struct MatchingControlStructure {
    232   bool foundflag;
    233   double bestL2;
    234   IndexList_t bestmatching;
    235   VectorArray_t oldpoints;
    236   VectorArray_t newpoints;
    237 };
    238 
    239 /** Recursive function to go through all possible matchings.
    240  *
    241  * \param _MCS structure holding global information to the recursion
    242  * \param _matching current matching being build up
    243  * \param _indices contains still available indices
    244  * \param _matchingsize
    245  */
    246 void recurseMatchings(
    247     MatchingControlStructure &_MCS,
    248     IndexList_t &_matching,
    249     IndexList_t _indices,
    250     unsigned int _matchingsize)
    251 {
    252   //!> threshold for L1 error below which matching is immediately acceptable
    253   const double L1THRESHOLD = 1e-2;
    254   if (!_MCS.foundflag) {
    255     if (_matching.size() < _matchingsize) {
    256       // go through all indices
    257       for (IndexList_t::iterator iter = _indices.begin();
    258           iter != _indices.end();) {
    259         // add index to matching
    260         _matching.push_back(*iter);
    261         // remove index but keep iterator to position (is the next to erase element)
    262         IndexList_t::iterator backupiter = _indices.erase(iter);
    263         // recurse with decreased _matchingsize
    264         recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
    265         // re-add chosen index and reset index to new position
    266         _indices.insert(backupiter, _matching.back());
    267         iter = backupiter;
    268         // remove index from _matching to make space for the next one
    269         _matching.pop_back();
    270       }
    271       // gone through all indices then exit recursion
    272       _MCS.foundflag = true;
    273     } else {
    274       // calculate errors
    275       std::pair<double, double> errors = calculateErrorOfMatching(
    276           _MCS.oldpoints, _MCS.newpoints, _matching);
    277       if (errors.first < L1THRESHOLD) {
    278         _MCS.bestmatching = _matching;
    279         _MCS.foundflag = true;
    280       }
    281       if (_MCS.bestL2 > errors.second) {
    282         _MCS.bestmatching = _matching;
    283         _MCS.bestL2 = errors.second;
    284       }
    285     }
    286   }
    287 }
    288 
    289 SphericalPointDistribution::Polygon_t
    290 SphericalPointDistribution::matchSphericalPointDistributions(
    291     const SphericalPointDistribution::Polygon_t &_polygon,
    292     const SphericalPointDistribution::Polygon_t &_newpolygon
    293     )
    294 {
    295   SphericalPointDistribution::Polygon_t remainingpoints;
    296   VectorArray_t remainingold(_polygon.begin(), _polygon.end());
    297   VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    298 
    299   if (_polygon.size() > 0) {
    300     MatchingControlStructure MCS;
    301     MCS.foundflag = false;
    302     MCS.bestL2 = std::numeric_limits<double>::max();
    303     MCS.oldpoints.insert(MCS.oldpoints.begin(), _polygon.begin(),_polygon.end() );
    304     MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
    305 
    306     // search for bestmatching combinatorially
    307     {
    308       // translate polygon into vector to enable index addressing
    309       IndexList_t indices(_newpolygon.size());
    310       std::generate(indices.begin(), indices.end(), UniqueNumber);
    311       IndexList_t matching;
    312 
    313       // walk through all matchings
    314       const unsigned int matchingsize = _polygon.size();
    315       ASSERT( matchingsize <= indices.size(),
    316           "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
    317       recurseMatchings(MCS, matching, indices, matchingsize);
    318     }
    319 
    320     // determine rotation angles to align the two point distributions with
    321     // respect to bestmatching
    322     std::vector<double> angles(3);
    323     Vector newCenter;
    324     {
    325       // calculate center of triangle/line/point consisting of first points of matching
    326       Vector oldCenter;
    327       IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    328       unsigned int i = 0;
    329       for (; (i<3) && (i<MCS.bestmatching.size()); ++i, ++iter) {
    330         oldCenter += remainingold[i];
    331         newCenter += remainingnew[*iter];
    332       }
    333       oldCenter *= 1./(double)i;
    334       newCenter *= 1./(double)i;
    335 
    336       Vector direction(0.,0.,0.);
    337       for(size_t i=0;i<NDIM;++i) {
    338         // create new rotation axis
    339         direction[i] = 1.;
    340         const Line axis (zeroVec, direction);
    341         // calculate rotation angle for this axis
    342         const double alpha = direction.Angle(oldCenter) - direction.Angle(newCenter);
    343         // perform rotation
    344         axis.rotateVector(newCenter, alpha);
    345         // store angle
    346         angles[i] = alpha;
    347         // reset direction component for next iteration
    348         direction[i] = 0.;
    349       }
    350     }
    351     const Line RotationAxis(zeroVec, newCenter);
    352     const double RotationAngle =
    353         newCenter.Angle(remainingold[0])
    354         - newCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
    355 
    356     // rotate _newpolygon
    357     SphericalPointDistribution::Polygon_t rotated_newpolygon =
    358         rotatePolygon(_newpolygon, angles, RotationAxis, RotationAngle);
    359 
    360     // remove all points in matching and return remaining ones
    361     return removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
    362   } else
    363     return _newpolygon;
    364 }
    365 
    366 template <>
    367 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>()
    368 {
    369   Polygon_t polygon;
    370   polygon.push_back( Vector(Bondlength,0.,0.));
    371   polygon.push_back( Vector(-Bondlength,0.,0.));
    372   ASSERT( polygon.size() == 2,
    373       "SphericalPointDistribution::get<2>() - polygon has wrong size.");
    374   return polygon;
    375 }
    376 
    377 template <>
    378 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>()
    379 {
    380   Polygon_t polygon;
    381   polygon.push_back( Vector(Bondlength,0.,0.));
    382   polygon.push_back( Vector(-Bondlength*0.5, SQRT_3*0.5,0.));
    383   polygon.push_back( Vector(-Bondlength*0.5, -SQRT_3*0.5,0.));
    384   ASSERT( polygon.size() == 3,
    385       "SphericalPointDistribution::get<3>() - polygon has wrong size.");
    386   return polygon;
    387 }
    388 
    389 template <>
    390 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>()
    391 {
    392   Polygon_t polygon;
    393   polygon.push_back( Vector(Bondlength,0.,0.));
    394   polygon.push_back( Vector(-Bondlength/3.0, Bondlength*2.0*M_SQRT2/3.0,0.));
    395   polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
    396   polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
    397   ASSERT( polygon.size() == 4,
    398       "SphericalPointDistribution::get<4>() - polygon has wrong size.");
    399   return polygon;
    400 }
    401 
    402 template <>
    403 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>()
    404 {
    405   Polygon_t polygon;
    406   polygon.push_back( Vector(Bondlength,0.,0.));
    407   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    408   polygon.push_back( Vector(0.0, Bondlength, 0.0));
    409   polygon.push_back( Vector(0.0, -Bondlength*0.5, Bondlength*SQRT_3*0.5));
    410   polygon.push_back( Vector(0.0, -Bondlength*0.5, -Bondlength*SQRT_3*0.5));
    411   ASSERT( polygon.size() == 5,
    412       "SphericalPointDistribution::get<5>() - polygon has wrong size.");
    413   return polygon;
    414 }
    415 
    416 template <>
    417 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>()
    418 {
    419   Polygon_t polygon;
    420   polygon.push_back( Vector(Bondlength,0.,0.));
    421   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    422   polygon.push_back( Vector(0.0, Bondlength, 0.0));
    423   polygon.push_back( Vector(0.0, -Bondlength, 0.0));
    424   polygon.push_back( Vector(0.0, 0.0, Bondlength));
    425   polygon.push_back( Vector(0.0, 0.0, -Bondlength));
    426   ASSERT( polygon.size() == 6,
    427       "SphericalPointDistribution::get<6>() - polygon has wrong size.");
    428   return polygon;
    429 }
    430 
    431 template <>
    432 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>()
    433 {
    434   Polygon_t polygon;
    435   polygon.push_back( Vector(Bondlength,0.,0.));
    436   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    437   polygon.push_back( Vector(0.0, Bondlength, 0.0));
    438   polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*0.4), Bondlength*sin(M_PI*0.4)));
    439   polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*0.8), Bondlength*sin(M_PI*0.8)));
    440   polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*1.2), Bondlength*sin(M_PI*1.2)));
    441   polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*1.6), Bondlength*sin(M_PI*1.6)));
    442   ASSERT( polygon.size() == 7,
    443       "SphericalPointDistribution::get<7>() - polygon has wrong size.");
    444   return polygon;
    445 }
    446 
    447 template <>
    448 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>()
    449 {
    450   Polygon_t polygon;
    451   polygon.push_back( Vector(Bondlength,0.,0.));
    452   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    453   polygon.push_back( Vector(-Bondlength/3.0, Bondlength*2.0*M_SQRT2/3.0, 0.0));
    454   polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
    455   polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
    456   polygon.push_back( Vector(Bondlength/3.0, -Bondlength*2.0*M_SQRT2/3.0, 0.0));
    457   polygon.push_back( Vector(Bondlength/3.0, Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
    458   polygon.push_back( Vector(Bondlength/3.0, Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
    459   ASSERT( polygon.size() == 8,
    460       "SphericalPointDistribution::get<8>() - polygon has wrong size.");
    461   return polygon;
    462 }
    463 
    464 template <>
    465 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>()
    466 {
    467   Polygon_t polygon;
    468   polygon.push_back( Vector(Bondlength,0.,0.));
    469   polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), Bondlength*sin(0.4*M_PI), 0.0));
    470   polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), -Bondlength*sin(0.4*M_PI), 0.0));
    471   polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), 0.0, Bondlength*sin(0.4*M_PI)));
    472   polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), 0.0, -Bondlength*sin(0.4*M_PI)));
    473   polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(0.25*M_PI), Bondlength*sin(0.8*M_PI)*cos(0.25*M_PI)));
    474   polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(0.75*M_PI), Bondlength*sin(0.8*M_PI)*cos(0.75*M_PI)));
    475   polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(1.25*M_PI), Bondlength*sin(0.8*M_PI)*cos(1.25*M_PI)));
    476   polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(1.75*M_PI), Bondlength*sin(0.8*M_PI)*cos(1.75*M_PI)));
    477   ASSERT( polygon.size() == 9,
    478       "SphericalPointDistribution::get<9>() - polygon has wrong size.");
    479   return polygon;
    480 }
    481 
    482 template <>
    483 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>()
    484 {
    485   Polygon_t polygon;
    486   polygon.push_back( Vector(Bondlength,0.,0.));
    487   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    488   const double temp_Bondlength = Bondlength*0.5*SQRT_3;
    489   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
    490   polygon.push_back( Vector(temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
    491   polygon.push_back( Vector(temp_Bondlength/SQRT_3, 0.0, temp_Bondlength));
    492   polygon.push_back( Vector(temp_Bondlength/SQRT_3, 0.0, -temp_Bondlength));
    493   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(0.25*M_PI), temp_Bondlength*cos(0.25*M_PI)));
    494   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(0.75*M_PI), temp_Bondlength*cos(0.75*M_PI)));
    495   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(1.25*M_PI), temp_Bondlength*cos(1.25*M_PI)));
    496   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(1.75*M_PI), temp_Bondlength*cos(1.75*M_PI)));
    497   ASSERT( polygon.size() == 10,
    498       "SphericalPointDistribution::get<10>() - polygon has wrong size.");
    499   return polygon;
    500 }
    501 
    502 template <>
    503 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>()
    504 {
    505   Polygon_t polygon;
    506   polygon.push_back( Vector(Bondlength,0.,0.));
    507   {
    508     const double temp_Bondlength = Bondlength*sin(0.4*M_PI);
    509     polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength, 0.0));
    510     polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(0.4*M_PI),  temp_Bondlength*sin(0.4*M_PI)));
    511     polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(0.8*M_PI),  temp_Bondlength*sin(0.8*M_PI)));
    512     polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(1.2*M_PI),  temp_Bondlength*sin(1.2*M_PI)));
    513     polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(1.6*M_PI),  temp_Bondlength*sin(1.6*M_PI)));
    514   }
    515   {
    516     const double temp_Bondlength = Bondlength*sin(0.8*M_PI);
    517     polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(0.2*M_PI),  temp_Bondlength*sin(0.2*M_PI)));
    518     polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(0.6*M_PI),  temp_Bondlength*sin(0.6*M_PI)));
    519     polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.0*M_PI),  temp_Bondlength*sin(1.0*M_PI)));
    520     polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.4*M_PI),  temp_Bondlength*sin(1.4*M_PI)));
    521     polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.8*M_PI),  temp_Bondlength*sin(1.8*M_PI)));
    522   }
    523   ASSERT( polygon.size() == 11,
    524       "SphericalPointDistribution::get<11>() - polygon has wrong size.");
    525   return polygon;
    526 }
    527 
    528 template <>
    529 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>()
    530 {
    531   Polygon_t polygon;
    532   polygon.push_back( Vector(Bondlength,0.,0.));
    533   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    534   const double temp_Bondlength = Bondlength*0.5*SQRT_3;
    535   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
    536   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(0.4*M_PI),  temp_Bondlength*sin(0.4*M_PI)));
    537   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(0.8*M_PI),  temp_Bondlength*sin(0.8*M_PI)));
    538   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(1.2*M_PI),  temp_Bondlength*sin(1.2*M_PI)));
    539   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(1.6*M_PI),  temp_Bondlength*sin(1.6*M_PI)));
    540   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
    541   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(0.4*M_PI),  -temp_Bondlength*sin(0.4*M_PI)));
    542   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(0.8*M_PI),  -temp_Bondlength*sin(0.8*M_PI)));
    543   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(1.2*M_PI),  -temp_Bondlength*sin(1.2*M_PI)));
    544   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(1.6*M_PI),  -temp_Bondlength*sin(1.6*M_PI)));
    545   ASSERT( polygon.size() == 12,
    546       "SphericalPointDistribution::get<12>() - polygon has wrong size.");
    547   return polygon;
    548 }
    549 
    550 template <>
    551 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>()
    552 {
    553   Polygon_t polygon;
    554   polygon.push_back( Vector(Bondlength,0.,0.));
    555   polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
    556   const double temp_Bondlength = Bondlength*0.5*SQRT_3;
    557   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
    558   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/3.0), temp_Bondlength*sin(M_PI/3.0)));
    559   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(2.0*M_PI/3.0), temp_Bondlength*sin(2.0*M_PI/3.0)));
    560   polygon.push_back( Vector(temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
    561   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(4.0*M_PI/3.0), temp_Bondlength*sin(4.0*M_PI/3.0)));
    562   polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(5.0*M_PI/3.0), temp_Bondlength*sin(5.0*M_PI/3.0)));
    563   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/6.0),  temp_Bondlength*sin(M_PI/6.0)));
    564   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(M_PI/3.0+M_PI/6.0)));
    565   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(2.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(2.0*M_PI/3.0 +M_PI/6.0)));
    566   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(3.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(3.0*M_PI/3.0 +M_PI/6.0)));
    567   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(4.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(4.0*M_PI/3.0 +M_PI/6.0)));
    568   polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(5.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(5.0*M_PI/3.0 +M_PI/6.0)));
    569   ASSERT( polygon.size() == 14,
    570       "SphericalPointDistribution::get<14>() - polygon has wrong size.");
    571   return polygon;
    572 }
     92template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>();
     93template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>();
     94template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>();
     95template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>();
     96template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>();
     97template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>();
     98template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>();
     99template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>();
     100template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>();
     101template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>();
     102template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>();
     103template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>();
     104template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>();
     105template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>();
    573106
    574107#endif /* SPHERICALPOINTDISTRIBUTION_HPP_ */

  • src/Fragmentation/Makefile.am

    r6ef48c r0d4daf  
    88        Fragmentation/Exporters/HydrogenPool.cpp \
    99        Fragmentation/Exporters/SaturatedFragment.cpp \
     10        Fragmentation/Exporters/SphericalPointDistribution.cpp \
     11        Fragmentation/Exporters/SphericalPointDistribution_getPoints.cpp \
    1012        Fragmentation/Homology/FragmentEdge.cpp \
    1113        Fragmentation/Homology/FragmentNode.cpp \
Note: See TracChangeset for help on using the changeset viewer.