Changeset 1bd79e


Ignore:
Timestamp:
Apr 15, 2010, 10:54:26 AM (15 years ago)
Author:
Tillmann Crueger <crueger@…>
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:
753f02
Parents:
273382
Message:

Changed implementation of Vector to forward operations to contained objects

Location:
src
Files:
21 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    r273382 r1bd79e  
    99                           gslvector.cpp \
    1010                           linearsystemofequations.cpp \
     11                           SingleVector.cpp \
    1112                           vector.cpp
    1213                           
    … …  
    1415                           gslvector.hpp \
    1516                           linearsystemofequations.hpp \
     17                           SingleVector.hpp \
    1618                           vector.hpp
    1719                           
    … …  
    122124                 periodentafel.cpp \
    123125                 Plane.cpp \
    124                  SingleVector.cpp \
    125126                 tesselation.cpp \
    126127                 tesselationhelpers.cpp \
    … …  
    158159         periodentafel.hpp \
    159160         Plane.hpp \
    160          SingleVector.hpp \
    161161         stackclass.hpp \
    162162         tesselation.hpp \

  • src/SingleVector.cpp

    r273382 r1bd79e  
    2121#include <gsl/gsl_vector.h>
    2222
     23#include <iostream>
     24
     25using namespace std;
     26
    2327/************************************ Functions for class vector ************************************/
    2428
    2529/** Constructor of class vector.
    2630 */
    27 SingleVector::SingleVector()
     31SingleVector::SingleVector() :
     32  Vector(Baseconstructor())
    2833{
    2934  x[0] = x[1] = x[2] = 0.;
    … …  
    3237/** Constructor of class vector.
    3338 */
    34 SingleVector::SingleVector(const double x1, const double x2, const double x3)
     39SingleVector::SingleVector(const double x1, const double x2, const double x3) :
     40  Vector(Baseconstructor())
    3541{
    3642  x[0] = x1;
    … …  
    4248 * Copy constructor
    4349 */
    44 SingleVector::SingleVector(const Vector& src)
     50SingleVector::SingleVector(const Vector& src) :
     51  Vector(Baseconstructor())
    4552{
    4653  x[0] = src[0];
    … …  
    4956}
    5057
    51 /**
    52  * Assignment operator
    53  */
    54 Vector& SingleVector::operator=(const Vector& src){
    55   // check for self assignment
    56   if(&src!=this){
    57     x[0] = src[0];
    58     x[1] = src[1];
    59     x[2] = src[2];
    60   }
    61   return *this;
     58SingleVector* SingleVector::clone() const{
     59  return new SingleVector(x[0],x[1],x[2]);
    6260}
    6361
    … …  
    7674    res += (x[i]-y[i])*(x[i]-y[i]);
    7775  return (res);
    78 };
    79 
    80 /** Calculates distance between this and another vector.
    81  * \param *y array to second vector
    82  * \return \f$| x - y |\f$
    83  */
    84 double SingleVector::Distance(const Vector &y) const
    85 {
    86   double res = 0.;
    87   for (int i=NDIM;i--;)
    88     res += (x[i]-y[i])*(x[i]-y[i]);
    89   return (sqrt(res));
    9076};
    9177
    … …  
    173159//  bool flag = false;
    174160  //vector Shifted, TranslationVector;
    175   Vector TestVector;
    176161//  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
    177162//  Log() << Verbose(2) << "Vector is: ";
    178163//  Output(out);
    179164//  Log() << Verbose(0) << endl;
    180   TestVector = *this;
     165  SingleVector TestVector(*this);
    181166  TestVector.InverseMatrixMultiplication(matrix);
    182167  for(int i=NDIM;i--;) { // correct periodically
    … …  
    208193
    209194
     195void SingleVector::AddVector(const Vector &y) {
     196  for(int i=NDIM;i--;)
     197    x[i] += y[i];
     198}
     199
     200void SingleVector::SubtractVector(const Vector &y){
     201  for(int i=NDIM;i--;)
     202    x[i] -= y[i];
     203}
     204
    210205/** Calculates VectorProduct between this and another vector.
    211206 *  -# returns the Product in place of vector from which it was initiated
    … …  
    216211void SingleVector::VectorProduct(const Vector &y)
    217212{
    218   Vector tmp;
     213  SingleVector tmp;
    219214  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
    220215  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
    221216  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
    222   (*this) = tmp;
     217  CopyVector(tmp);
    223218};
    224219
    … …  
    245240double SingleVector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
    246241{
    247   Vector temp;
    248 
    249242  // first create part that is orthonormal to PlaneNormal with withdraw
    250   temp = (*this )- PlaneOffset;
     243  Vector temp = VecFromRep(this)- PlaneOffset;
    251244  temp.MakeNormalTo(PlaneNormal);
    252245  temp.Scale(-1.);
    253246  // then add connecting vector from plane to point
    254   temp += (*this)-PlaneOffset;
     247  temp += VecFromRep(this)-PlaneOffset;
    255248  double sign = temp.ScalarProduct(PlaneNormal);
    256249  if (fabs(sign) > MYEPSILON)
    … …  
    282275
    283276  return helper;
    284 };
    285 
    286 /** Calculates norm of this vector.
    287  * \return \f$|x|\f$
    288  */
    289 double SingleVector::Norm() const
    290 {
    291   double res = 0.;
    292   for (int i=NDIM;i--;)
    293     res += this->x[i]*this->x[i];
    294   return (sqrt(res));
    295 };
    296 
    297 /** Calculates squared norm of this vector.
    298  * \return \f$|x|^2\f$
    299  */
    300 double SingleVector::NormSquared() const
    301 {
    302   return (ScalarProduct(*this));
    303 };
    304 
    305 /** Normalizes this vector.
    306  */
    307 void SingleVector::Normalize()
    308 {
    309   double res = 0.;
    310   for (int i=NDIM;i--;)
    311     res += this->x[i]*this->x[i];
    312   if (fabs(res) > MYEPSILON)
    313     res = 1./sqrt(res);
    314   Scale(&res);
    315 };
    316 
    317 /** Zeros all components of this vector.
    318  */
    319 void SingleVector::Zero()
    320 {
    321   for (int i=NDIM;i--;)
    322     this->x[i] = 0.;
    323 };
    324 
    325 /** Zeros all components of this vector.
    326  */
    327 void SingleVector::One(const double one)
    328 {
    329   for (int i=NDIM;i--;)
    330     this->x[i] = one;
    331 };
    332 
    333 /** Initialises all components of this vector.
    334  */
    335 void SingleVector::Init(const double x1, const double x2, const double x3)
    336 {
    337   x[0] = x1;
    338   x[1] = x2;
    339   x[2] = x3;
    340277};
    341278
    … …  
    410347}
    411348
    412 double& SingleVector::at(size_t i){
    413   return (*this)[i];
    414 }
    415 
    416 const double& SingleVector::at(size_t i) const{
    417   return (*this)[i];
    418 }
    419 
    420349double* SingleVector::get(){
    421350  return x;
    … …  
    424353 * \param *factor pointer to scaling factor
    425354 */
    426 void SingleVector::Scale(const double ** const factor)
    427 {
    428   for (int i=NDIM;i--;)
    429     x[i] *= (*factor)[i];
    430 };
    431 
    432 void SingleVector::Scale(const double * const factor)
    433 {
    434   for (int i=NDIM;i--;)
    435     x[i] *= *factor;
     355void SingleVector::ScaleAll(const double *factor)
     356{
     357  for (int i=NDIM;i--;)
     358    x[i] *= factor[i];
    436359};
    437360
    … …  
    440363  for (int i=NDIM;i--;)
    441364    x[i] *= factor;
    442 };
    443 
    444 /** Translate atom by given vector.
    445  * \param trans[] translation vector.
    446  */
    447 void SingleVector::Translate(const Vector &trans)
    448 {
    449   for (int i=NDIM;i--;)
    450     x[i] += trans[i];
    451365};
    452366
    … …  
    474388void SingleVector::MatrixMultiplication(const double * const M)
    475389{
    476   Vector C;
     390  SingleVector C;
    477391  // do the matrix multiplication
    478392  C[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    … …  
    480394  C[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    481395
    482   *this = C;
     396  CopyVector(C);
    483397};
    484398
    … …  
    488402bool SingleVector::InverseMatrixMultiplication(const double * const A)
    489403{
    490   Vector C;
     404  SingleVector C;
    491405  double B[NDIM*NDIM];
    492406  double detA = RDET3(A);
    … …  
    511425    C[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    512426
    513     *this = C;
     427    CopyVector(C);
    514428    return true;
    515429  } else {
    … …  
    518432};
    519433
    520 
    521 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
    522  * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
    523  * \param *x1 first vector
    524  * \param *x2 second vector
    525  * \param *x3 third vector
    526  * \param *factors three-component vector with the factor for each given vector
    527  */
    528 void SingleVector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
    529 {
    530   for(int i=NDIM;i--;)
    531     x[i] = factors[0]*x1[i] + factors[1]*x2[i] + factors[2]*x3[i];
    532 };
    533434
    534435/** Mirrors atom against a given plane.
    … …  
    615516bool SingleVector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
    616517{
    617   Vector a = (*this) - offset;
     518  Vector a = VecFromRep(this);
     519  a -= offset;
    618520  a.InverseMatrixMultiplication(parallelepiped);
    619521  bool isInside = true;
    … …  
    624526  return isInside;
    625527}
     528
     529
     530void SingleVector::CopyVector(SingleVector& vec) {
     531  for(int i =NDIM; i--;)
     532    x[i]=vec.x[i];
     533}
     534
     535bool SingleVector::isBaseClass() const{
     536  return false;
     537}

  • src/SingleVector.hpp

    r273382 r1bd79e  
    3131  virtual ~SingleVector();
    3232
    33   virtual double Distance(const Vector &y) const;
    3433  virtual double DistanceSquared(const Vector &y) const;
    3534  virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
    … …  
    3736  virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
    3837  virtual double ScalarProduct(const Vector &y) const;
    39   virtual double Norm() const;
    40   virtual double NormSquared() const;
    4138  virtual double Angle(const Vector &y) const;
    4239  virtual bool IsZero() const;
    … …  
    4542  virtual bool IsEqualTo(const Vector &a) const;
    4643
     44  virtual void AddVector(const Vector &y);
     45  virtual void SubtractVector(const Vector &y);
    4746  virtual void VectorProduct(const Vector &y);
    4847  virtual void ProjectOntoPlane(const Vector &y);
    4948  virtual void ProjectIt(const Vector &y);
    5049  virtual Vector Projection(const Vector &y) const;
    51   virtual void Zero();
    52   virtual void One(const double one);
    53   virtual void Init(const double x1, const double x2, const double x3);
    54   virtual void Normalize();
    55   virtual void Translate(const Vector &x);
    5650  virtual void Mirror(const Vector &x);
    57   virtual void Scale(const double ** const factor);
    58   virtual void Scale(const double * const factor);
     51  virtual void ScaleAll(const double *factor);
    5952  virtual void Scale(const double factor);
    6053  virtual void MatrixMultiplication(const double * const M);
    6154  virtual bool InverseMatrixMultiplication(const double * const M);
    6255  virtual void KeepPeriodic(const double * const matrix);
    63   virtual void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
    6456  virtual bool GetOneNormalVector(const Vector &x1);
    6557  virtual bool MakeNormalTo(const Vector &y1);
    66   //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
    6758  virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
    6859  virtual void WrapPeriodically(const double * const M, const double * const Minv);
    … …  
    7162  virtual double& operator[](size_t i);
    7263  virtual const double& operator[](size_t i) const;
    73   virtual double& at(size_t i);
    74   virtual const double& at(size_t i) const;
    75 
    76   // Assignment operator
    77   virtual Vector &operator=(const Vector& src);
    7864
    7965  // Access to internal structure
    8066  virtual double* get();
     67protected:
    8168
    8269private:
     70  // method used for protection, i.e. to avoid infinite recursion
     71  // when our internal rep becomes messed up
     72  virtual bool isBaseClass() const;
     73  virtual SingleVector *clone() const;
     74  void CopyVector(SingleVector &rhs);
    8375  double x[NDIM];
    8476

  • src/boundary.cpp

    r273382 r1bd79e  
    827827
    828828  // calculate filler grid in [0,1]^3
    829   FillerDistance.Init(distance[0], distance[1], distance[2]);
     829  FillerDistance = Vector(distance[0], distance[1], distance[2]);
    830830  FillerDistance.InverseMatrixMultiplication(M);
    831831  for(int i=0;i<NDIM;i++)
    … …  
    841841      for (n[2] = 0; n[2] < N[2]; n[2]++) {
    842842        // calculate position of current grid vector in untransformed box
    843         CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     843        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    844844        CurrentPosition.MatrixMultiplication(M);
    845845        Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;

  • src/builder.cpp

    r273382 r1bd79e  
    16441644                first = World::getInstance().createAtom();
    16451645                first->type = periode->FindElement(1);
    1646                 first->x.Init(0.441, -0.143, 0.);
     1646                first->x = Vector(0.441, -0.143, 0.);
    16471647                filler->AddAtom(first);
    16481648                second = World::getInstance().createAtom();
    16491649                second->type = periode->FindElement(1);
    1650                 second->x.Init(-0.464, 1.137, 0.0);
     1650                second->x = Vector(-0.464, 1.137, 0.0);
    16511651                filler->AddAtom(second);
    16521652                third = World::getInstance().createAtom();
    16531653                third->type = periode->FindElement(8);
    1654                 third->x.Init(-0.464, 0.177, 0.);
     1654                third->x = Vector(-0.464, 0.177, 0.);
    16551655                filler->AddAtom(third);
    16561656                filler->AddBond(first, third, 1);

  • src/ellipsoid.cpp

    r273382 r1bd79e  
    5959  Matrix[8] = cos(theta);
    6060  helper.MatrixMultiplication(Matrix);
    61   helper.Scale(InverseLength);
     61  helper.ScaleAll(InverseLength);
    6262  //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
    6363
    … …  
    7070  theta = -EllipsoidAngle[1];
    7171  phi = -EllipsoidAngle[2];
    72   helper.Scale(EllipsoidLength);
     72  helper.ScaleAll(EllipsoidLength);
    7373  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
    7474  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);

  • src/molecule.cpp

    r273382 r1bd79e  
    239239    matrix = ReturnFullMatrixforSymmetric(cell_size);
    240240    Orthovector1.MatrixMultiplication(matrix);
    241     InBondvector.SubtractVector(Orthovector1); // subtract just the additional translation
     241    InBondvector -= Orthovector1; // subtract just the additional translation
    242242    Free(&matrix);
    243243    bondlength = InBondvector.Norm();
    … …  
    272272        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
    273273      }
    274       InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
     274      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
    275275      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
    276276      FirstOtherAtom->x = InBondvector;  // ... and add distance vector to replacement atom
    … …  
    356356        SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
    357357      }
    358       FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
    359       SecondOtherAtom->x.Scale(&BondRescale);
     358      FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
     359      SecondOtherAtom->x *= BondRescale;
    360360      //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
    361361      for(int i=NDIM;i--;) { // and make relative to origin atom

  • src/molecule_geometry.cpp

    r273382 r1bd79e  
    191191/** Scales all atoms by \a *factor.
    192192 * \param *factor pointer to scaling factor
     193 *
     194 * TODO: Is this realy what is meant, i.e.
     195 * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
     196 * or rather
     197 * x=(**factor) * x (as suggested by comment)
    193198 */
    194199void molecule::Scale(const double ** const factor)
    … …  
    199204    ptr = ptr->next;
    200205    for (int j=0;j<MDSteps;j++)
    201       ptr->Trajectory.R.at(j).Scale(factor);
    202     ptr->x.Scale(factor);
     206      ptr->Trajectory.R.at(j).ScaleAll(*factor);
     207    ptr->x.ScaleAll(*factor);
    203208  }
    204209};
    … …  
    214219    ptr = ptr->next;
    215220    for (int j=0;j<MDSteps;j++)
    216       ptr->Trajectory.R.at(j).Translate(*trans);
    217     ptr->x.Translate(*trans);
     221      ptr->Trajectory.R.at(j) += (*trans);
     222    ptr->x += (*trans);
    218223  }
    219224};

  • src/tesselationhelpers.cpp

    r273382 r1bd79e  
    427427  double volume;
    428428
    429   TetraederVector[0].CopyVector(a);
    430   TetraederVector[1].CopyVector(b);
    431   TetraederVector[2].CopyVector(c);
     429  TetraederVector[0] = a;
     430  TetraederVector[1] = b;
     431  TetraederVector[2] = c;
    432432  for (int j=0;j<3;j++)
    433433    TetraederVector[j].SubtractVector(d);
    434   Point.CopyVector(TetraederVector[0]);
     434  Point = TetraederVector[0];
    435435  Point.VectorProduct(TetraederVector[1]);
    436436  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));

  • src/unittests/ActOnAllUnitTest.cpp

    r273382 r1bd79e  
    7070
    7171  // scaling by value
    72   VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 2. );
     72  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, factor );
    7373  CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
    7474
    75   VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 0.5 );
    76   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
    77 
    78   // scaling by ref
    79   VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&factor );
    80   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
    81 
    82   VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&inverse );
     75  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, inverse );
    8376  CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
    8477
    … …  
    9083    inverses[i] = 1./factors[i];
    9184  }
    92   VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&factors );
     85  VL.ActOnAll<Vector,void,const double*>(&Vector::ScaleAll, factors );
    9386  CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
    9487
    95   VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&inverses );
     88  VL.ActOnAll<Vector,void,const double*>(&Vector::ScaleAll, inverses );
    9689  CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
    9790  Free(factors);

  • src/unittests/AnalysisCorrelationToPointUnitTest.cpp

    r273382 r1bd79e  
    6464  Walker = World::getInstance().createAtom();
    6565  Walker->type = hydrogen;
    66   Walker->node->Init(1., 0., 1. );
     66  *Walker->node = Vector(1., 0., 1. );
    6767  TestMolecule->AddAtom(Walker);
    6868  Walker = World::getInstance().createAtom();
    6969  Walker->type = hydrogen;
    70   Walker->node->Init(0., 1., 1. );
     70  *Walker->node = Vector(0., 1., 1. );
    7171  TestMolecule->AddAtom(Walker);
    7272  Walker = World::getInstance().createAtom();
    7373  Walker->type = hydrogen;
    74   Walker->node->Init(1., 1., 0. );
     74  *Walker->node = Vector(1., 1., 0. );
    7575  TestMolecule->AddAtom(Walker);
    7676  Walker = World::getInstance().createAtom();
    7777  Walker->type = hydrogen;
    78   Walker->node->Init(0., 0., 0. );
     78  *Walker->node = Vector(0., 0., 0. );
    7979  TestMolecule->AddAtom(Walker);
    8080

  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

    r273382 r1bd79e  
    4040void AnalysisCorrelationToSurfaceUnitTest::setUp()
    4141{
    42   ASSERT_DO(Assert::Throw);
     42  //ASSERT_DO(Assert::Throw);
    4343
    4444  atom *Walker = NULL;
    … …  
    7373  Walker = World::getInstance().createAtom();
    7474  Walker->type = hydrogen;
    75   Walker->node->Init(1., 0., 1. );
    76   TestMolecule->AddAtom(Walker);
    77   Walker = World::getInstance().createAtom();
    78   Walker->type = hydrogen;
    79   Walker->node->Init(0., 1., 1. );
    80   TestMolecule->AddAtom(Walker);
    81   Walker = World::getInstance().createAtom();
    82   Walker->type = hydrogen;
    83   Walker->node->Init(1., 1., 0. );
    84   TestMolecule->AddAtom(Walker);
    85   Walker = World::getInstance().createAtom();
    86   Walker->type = hydrogen;
    87   Walker->node->Init(0., 0., 0. );
     75  *Walker->node = Vector(1., 0., 1. );
     76  TestMolecule->AddAtom(Walker);
     77  Walker = World::getInstance().createAtom();
     78  Walker->type = hydrogen;
     79  *Walker->node = Vector(0., 1., 1. );
     80  TestMolecule->AddAtom(Walker);
     81  Walker = World::getInstance().createAtom();
     82  Walker->type = hydrogen;
     83  *Walker->node = Vector(1., 1., 0. );
     84  TestMolecule->AddAtom(Walker);
     85  Walker = World::getInstance().createAtom();
     86  Walker->type = hydrogen;
     87  *Walker->node = Vector(0., 0., 0. );
    8888  TestMolecule->AddAtom(Walker);
    8989
    … …  
    106106  Walker = World::getInstance().createAtom();
    107107  Walker->type = carbon;
    108   Walker->node->Init(4., 0., 4. );
    109   TestMolecule->AddAtom(Walker);
    110   Walker = World::getInstance().createAtom();
    111   Walker->type = carbon;
    112   Walker->node->Init(0., 4., 4. );
    113   TestMolecule->AddAtom(Walker);
    114   Walker = World::getInstance().createAtom();
    115   Walker->type = carbon;
    116   Walker->node->Init(4., 4., 0. );
     108  *Walker->node = Vector(4., 0., 4. );
     109  TestMolecule->AddAtom(Walker);
     110  Walker = World::getInstance().createAtom();
     111  Walker->type = carbon;
     112  *Walker->node = Vector(0., 4., 4. );
     113  TestMolecule->AddAtom(Walker);
     114  Walker = World::getInstance().createAtom();
     115  Walker->type = carbon;
     116  *Walker->node = Vector(4., 4., 0. );
    117117  TestMolecule->AddAtom(Walker);
    118118  // add inner atoms
    119119  Walker = World::getInstance().createAtom();
    120120  Walker->type = carbon;
    121   Walker->node->Init(0.5, 0.5, 0.5 );
     121  *Walker->node = Vector(0.5, 0.5, 0.5 );
    122122  TestMolecule->AddAtom(Walker);
    123123

  • src/unittests/AnalysisPairCorrelationUnitTest.cpp

    r273382 r1bd79e  
    6363  Walker = World::getInstance().createAtom();
    6464  Walker->type = hydrogen;
    65   Walker->node->Init(1., 0., 1. );
     65  *Walker->node = Vector(1., 0., 1. );
    6666  TestMolecule->AddAtom(Walker);
    6767  Walker = World::getInstance().createAtom();
    6868  Walker->type = hydrogen;
    69   Walker->node->Init(0., 1., 1. );
     69  *Walker->node = Vector(0., 1., 1. );
    7070  TestMolecule->AddAtom(Walker);
    7171  Walker = World::getInstance().createAtom();
    7272  Walker->type = hydrogen;
    73   Walker->node->Init(1., 1., 0. );
     73  *Walker->node = Vector(1., 1., 0. );
    7474  TestMolecule->AddAtom(Walker);
    7575  Walker = World::getInstance().createAtom();
    7676  Walker->type = hydrogen;
    77   Walker->node->Init(0., 0., 0. );
     77  *Walker->node = Vector(0., 0., 0. );
    7878  TestMolecule->AddAtom(Walker);
    7979

  • src/unittests/analysisbondsunittest.cpp

    r273382 r1bd79e  
    6969  Walker = World::getInstance().createAtom();
    7070  Walker->type = hydrogen;
    71   Walker->node->Init(1.5, 0., 1.5 );
     71  *Walker->node = Vector(1.5, 0., 1.5 );
    7272  TestMolecule->AddAtom(Walker);
    7373  Walker = World::getInstance().createAtom();
    7474  Walker->type = hydrogen;
    75   Walker->node->Init(0., 1.5, 1.5 );
     75  *Walker->node = Vector(0., 1.5, 1.5 );
    7676  TestMolecule->AddAtom(Walker);
    7777  Walker = World::getInstance().createAtom();
    7878  Walker->type = hydrogen;
    79   Walker->node->Init(1.5, 1.5, 0. );
     79  *Walker->node = Vector(1.5, 1.5, 0. );
    8080  TestMolecule->AddAtom(Walker);
    8181  Walker = World::getInstance().createAtom();
    8282  Walker->type = hydrogen;
    83   Walker->node->Init(0., 0., 0. );
     83  *Walker->node = Vector(0., 0., 0. );
    8484  TestMolecule->AddAtom(Walker);
    8585  Walker = World::getInstance().createAtom();
    8686  Walker->type = carbon;
    87   Walker->node->Init(0.5, 0.5, 0.5 );
     87  *Walker->node = Vector(0.5, 0.5, 0.5 );
    8888  TestMolecule->AddAtom(Walker);
    8989

  • src/unittests/bondgraphunittest.cpp

    r273382 r1bd79e  
    6565  Walker = World::getInstance().createAtom();
    6666  Walker->type = hydrogen;
    67   Walker->node->Init(1., 0., 1. );
     67  *Walker->node = Vector(1., 0., 1. );
    6868  TestMolecule->AddAtom(Walker);
    6969  Walker = World::getInstance().createAtom();
    7070  Walker->type = hydrogen;
    71   Walker->node->Init(0., 1., 1. );
     71  *Walker->node = Vector(0., 1., 1. );
    7272  TestMolecule->AddAtom(Walker);
    7373  Walker = World::getInstance().createAtom();
    7474  Walker->type = hydrogen;
    75   Walker->node->Init(1., 1., 0. );
     75  *Walker->node = Vector(1., 1., 0. );
    7676  TestMolecule->AddAtom(Walker);
    7777  Walker = World::getInstance().createAtom();
    7878  Walker->type = hydrogen;
    79   Walker->node->Init(0., 0., 0. );
     79  *Walker->node = Vector(0., 0., 0. );
    8080  TestMolecule->AddAtom(Walker);
    8181

  • src/unittests/listofbondsunittest.cpp

    r273382 r1bd79e  
    5858  Walker = World::getInstance().createAtom();
    5959  Walker->type = hydrogen;
    60   Walker->node->Init(1., 0., 1. );
    61   TestMolecule->AddAtom(Walker);
    62   Walker = World::getInstance().createAtom();
    63   Walker->type = hydrogen;
    64   Walker->node->Init(0., 1., 1. );
    65   TestMolecule->AddAtom(Walker);
    66   Walker = World::getInstance().createAtom();
    67   Walker->type = hydrogen;
    68   Walker->node->Init(1., 1., 0. );
    69   TestMolecule->AddAtom(Walker);
    70   Walker = World::getInstance().createAtom();
    71   Walker->type = hydrogen;
    72   Walker->node->Init(0., 0., 0. );
     60  *Walker->node = Vector(1., 0., 1. );
     61  TestMolecule->AddAtom(Walker);
     62  Walker = World::getInstance().createAtom();
     63  Walker->type = hydrogen;
     64  *Walker->node = Vector(0., 1., 1. );
     65  TestMolecule->AddAtom(Walker);
     66  Walker = World::getInstance().createAtom();
     67  Walker->type = hydrogen;
     68  *Walker->node = Vector(1., 1., 0. );
     69  TestMolecule->AddAtom(Walker);
     70  Walker = World::getInstance().createAtom();
     71  Walker->type = hydrogen;
     72  *Walker->node = Vector(0., 0., 0. );
    7373  TestMolecule->AddAtom(Walker);
    7474

  • src/unittests/tesselation_boundarytriangleunittest.cpp

    r273382 r1bd79e  
    8686
    8787  // simple test on y line
    88   Point.Init(-1.,0.5,0.);
    89   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    90   Point.Init(0.,0.5,0.);
    91   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    92   Point.Init(-4.,0.5,0.);
     88  Point = Vector(-1.,0.5,0.);
     89  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
     90  Point = Vector(0.,0.5,0.);
     91  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     92  Point = Vector(-4.,0.5,0.);
    9393  CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    94   Point.Init(0.,0.5,0.);
     94  Point = Vector(0.,0.5,0.);
    9595  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    9696
    9797  // simple test on x line
    98   Point.Init(0.5,-1.,0.);
    99   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    100   Point.Init(0.5,0.,0.);
    101   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    102   Point.Init(0.5,-6.,0.);
     98  Point = Vector(0.5,-1.,0.);
     99  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
     100  Point = Vector(0.5,0.,0.);
     101  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     102  Point = Vector(0.5,-6.,0.);
    103103  CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    104   Point.Init(0.5,0.,0.);
     104  Point = Vector(0.5,0.,0.);
    105105  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    106106
    107107  // simple test on slanted line
    108   Point.Init(1.,1.,0.);
     108  Point = Vector(1.,1.,0.);
    109109  CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    110   Point.Init(0.5,0.5,0.);
    111   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    112   Point.Init(5.,5.,0.);
     110  Point = Vector(0.5,0.5,0.);
     111  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     112  Point = Vector(5.,5.,0.);
    113113  CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    114   Point.Init(0.5,0.5,0.);
     114  Point = Vector(0.5,0.5,0.);
    115115  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    116116
    117117  // simple test on first node
    118   Point.Init(-1.,-1.,0.);
     118  Point = Vector(-1.,-1.,0.);
    119119  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    120   Point.Init(0.,0.,0.);
     120  Point = Vector(0.,0.,0.);
    121121  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    122122
    123123  // simple test on second node
    124   Point.Init(0.,2.,0.);
    125   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    126   Point.Init(0.,1.,0.);
     124  Point = Vector(0.,2.,0.);
     125  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
     126  Point = Vector(0.,1.,0.);
    127127  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    128128
    129129  // simple test on third node
    130   Point.Init(2.,0.,0.);
    131   CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    132   Point.Init(1.,0.,0.);
     130  Point = Vector(2.,0.,0.);
     131  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
     132  Point = Vector(1.,0.,0.);
    133133  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    134134};
    … …  
    142142
    143143  // straight down/up
    144   Point.Init(1./3.,1./3.,+5.);
     144  Point = Vector(1./3.,1./3.,+5.);
    145145  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    146   Point.Init(1./3.,1./3.,0.);
    147   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    148   Point.Init(1./3.,1./3.,-5.);
     146  Point = Vector(1./3.,1./3.,0.);
     147  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     148  Point = Vector(1./3.,1./3.,-5.);
    149149  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    150   Point.Init(1./3.,1./3.,0.);
     150  Point = Vector(1./3.,1./3.,0.);
    151151  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    152152
    153153  // simple test on y line
    154   Point.Init(-1.,0.5,+2.);
     154  Point = Vector(-1.,0.5,+2.);
    155155  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    156   Point.Init(0.,0.5,0.);
    157   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    158   Point.Init(-1.,0.5,-3.);
     156  Point = Vector(0.,0.5,0.);
     157  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     158  Point = Vector(-1.,0.5,-3.);
    159159  CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    160   Point.Init(0.,0.5,0.);
     160  Point = Vector(0.,0.5,0.);
    161161  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    162162
    163163  // simple test on x line
    164   Point.Init(0.5,-1.,+1.);
     164  Point = Vector(0.5,-1.,+1.);
    165165  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    166   Point.Init(0.5,0.,0.);
    167   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    168   Point.Init(0.5,-1.,-2.);
     166  Point = Vector(0.5,0.,0.);
     167  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     168  Point = Vector(0.5,-1.,-2.);
    169169  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    170   Point.Init(0.5,0.,0.);
     170  Point = Vector(0.5,0.,0.);
    171171  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    172172
    173173  // simple test on slanted line
    174   Point.Init(1.,1.,+3.);
     174  Point = Vector(1.,1.,+3.);
    175175  CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    176   Point.Init(0.5,0.5,0.);
    177   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    178   Point.Init(1.,1.,-4.);
     176  Point = Vector(0.5,0.5,0.);
     177  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     178  Point = Vector(1.,1.,-4.);
    179179  CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    180   Point.Init(0.5,0.5,0.);
     180  Point = Vector(0.5,0.5,0.);
    181181  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    182182
    183183  // simple test on first node
    184   Point.Init(-1.,-1.,5.);
     184  Point = Vector(-1.,-1.,5.);
    185185  CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    186   Point.Init(0.,0.,0.);
     186  Point = Vector(0.,0.,0.);
    187187  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    188188
    189189  // simple test on second node
    190   Point.Init(0.,2.,5.);
     190  Point = Vector(0.,2.,5.);
    191191  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    192   Point.Init(0.,1.,0.);
     192  Point = Vector(0.,1.,0.);
    193193  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    194194
    195195  // simple test on third node
    196   Point.Init(2.,0.,5.);
     196  Point = Vector(2.,0.,5.);
    197197  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
    198   Point.Init(1.,0.,0.);
    199   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
    200 };
     198  Point = Vector(1.,0.,0.);
     199  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
     200};

  • src/unittests/vectorunittest.cpp

    r273382 r1bd79e  
    2626#endif /*HAVE_TESTRUNNER*/
    2727
     28#include <iostream>
     29
     30using namespace std;
     31
    2832/********************************************** Test classes **************************************/
    2933
    … …  
    3438void VectorTest::setUp()
    3539{
    36   zero.Init(0.,0.,0.);
    37   unit.Init(1.,0.,0.);
    38   otherunit.Init(0.,1.,0.);
    39   notunit.Init(0.,1.,1.);
    40   two.Init(2.,1.,0.);
     40  zero  = Vector(0.,0.,0.);
     41  unit = Vector(1.,0.,0.);
     42  otherunit = Vector(0.,1.,0.);
     43  notunit = Vector(0.,1.,1.);
     44  two = Vector(2.,1.,0.);
    4145};
    4246
    … …  
    6973  double factor;
    7074  // copy vector
    71   fixture.Init(2.,3.,4.);
     75  fixture = Vector(2.,3.,4.);
    7276  CPPUNIT_ASSERT_EQUAL( Vector(2.,3.,4.), fixture );
    7377  // summation and scaling
    … …  
    223227void VectorTest::VectorRotationTest()
    224228{
    225   fixture.Init(-1.,0.,0.);
     229  fixture = Vector(-1.,0.,0.);
    226230
    227231  // zero vector does not change
    … …  
    266270  parallelepiped[8] = 1;
    267271
    268   fixture.CopyVector(zero);
    269   CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    270   fixture.Init(2.5,2.5,2.5);
     272  fixture = zero;
     273  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     274  fixture = Vector(2.5,2.5,2.5);
    271275  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    272   fixture.Init(1.,1.,1.);
    273   CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    274   fixture.Init(3.5,3.5,3.5);
    275   CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    276   fixture.Init(2.,2.,2.);
     276  fixture = Vector(1.,1.,1.);
     277  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     278  fixture = Vector(3.5,3.5,3.5);
     279  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     280  fixture = Vector(2.,2.,2.);
    277281  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    278   fixture.Init(2.,3.,2.);
     282  fixture = Vector(2.,3.,2.);
    279283  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    280   fixture.Init(-2.,2.,-1.);
    281   CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
    282 }
    283 
     284  fixture = Vector(-2.,2.,-1.);
     285  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     286}
     287

  • src/vector.cpp

    r273382 r1bd79e  
    66
    77
    8 #include "defs.hpp"
    9 #include "gslmatrix.hpp"
    10 #include "leastsquaremin.hpp"
    11 #include "memoryallocator.hpp"
    12 #include "vector.hpp"
    13 #include "Helpers/fast_functions.hpp"
     8#include "SingleVector.hpp"
    149#include "Helpers/Assert.hpp"
    15 #include "Plane.hpp"
    16 #include "Exceptions/LinearDependenceException.hpp"
    17 
    18 #include <gsl/gsl_linalg.h>
    19 #include <gsl/gsl_matrix.h>
    20 #include <gsl/gsl_permutation.h>
    21 #include <gsl/gsl_vector.h>
     10
     11#include <iostream>
     12
     13using namespace std;
     14
    2215
    2316/************************************ Functions for class vector ************************************/
    … …  
    2518/** Constructor of class vector.
    2619 */
    27 Vector::Vector()
    28 {
    29   x[0] = x[1] = x[2] = 0.;
    30 };
     20Vector::Vector() :
     21  rep(new SingleVector())
     22{};
     23
     24Vector::Vector(Baseconstructor)  // used by derived objects to construct their bases
     25{}
     26
     27Vector::Vector(Baseconstructor,const Vector* v) :
     28  rep(v->clone())
     29{}
     30
     31Vector Vector::VecFromRep(const Vector* v){
     32  return Vector(Baseconstructor(),v);
     33}
    3134
    3235/** Constructor of class vector.
    3336 */
    34 Vector::Vector(const double x1, const double x2, const double x3)
    35 {
    36   x[0] = x1;
    37   x[1] = x2;
    38   x[2] = x3;
    39 };
     37Vector::Vector(const double x1, const double x2, const double x3) :
     38  rep(new SingleVector(x1,x2,x3))
     39{};
    4040
    4141/**
    4242 * Copy constructor
    4343 */
    44 Vector::Vector(const Vector& src)
    45 {
    46   x[0] = src[0];
    47   x[1] = src[1];
    48   x[2] = src[2];
    49 }
     44Vector::Vector(const Vector& src) :
     45  rep(src.rep->clone())
     46{}
    5047
    5148/**
    … …  
    5350 */
    5451Vector& Vector::operator=(const Vector& src){
     52  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
    5553  // check for self assignment
    5654  if(&src!=this){
    57     x[0] = src[0];
    58     x[1] = src[1];
    59     x[2] = src[2];
     55    rep.reset(src.rep->clone());
    6056  }
    6157  return *this;
    … …  
    7268double Vector::DistanceSquared(const Vector &y) const
    7369{
    74   double res = 0.;
    75   for (int i=NDIM;i--;)
    76     res += (x[i]-y[i])*(x[i]-y[i]);
    77   return (res);
     70  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     71  return rep->DistanceSquared(y);
    7872};
    7973
    … …  
    9488double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
    9589{
    96   double res = Distance(y), tmp, matrix[NDIM*NDIM];
    97   Vector Shiftedy, TranslationVector;
    98   int N[NDIM];
    99   matrix[0] = cell_size[0];
    100   matrix[1] = cell_size[1];
    101   matrix[2] = cell_size[3];
    102   matrix[3] = cell_size[1];
    103   matrix[4] = cell_size[2];
    104   matrix[5] = cell_size[4];
    105   matrix[6] = cell_size[3];
    106   matrix[7] = cell_size[4];
    107   matrix[8] = cell_size[5];
    108   // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    109   for (N[0]=-1;N[0]<=1;N[0]++)
    110     for (N[1]=-1;N[1]<=1;N[1]++)
    111       for (N[2]=-1;N[2]<=1;N[2]++) {
    112         // create the translation vector
    113         TranslationVector.Zero();
    114         for (int i=NDIM;i--;)
    115           TranslationVector.x[i] = (double)N[i];
    116         TranslationVector.MatrixMultiplication(matrix);
    117         // add onto the original vector to compare with
    118         Shiftedy = y + TranslationVector;
    119         // get distance and compare with minimum so far
    120         tmp = Distance(Shiftedy);
    121         if (tmp < res) res = tmp;
    122       }
    123   return (res);
     90  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     91  return rep->PeriodicDistance(y,cell_size);
    12492};
    12593
    … …  
    13199double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
    132100{
    133   double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    134   Vector Shiftedy, TranslationVector;
    135   int N[NDIM];
    136   matrix[0] = cell_size[0];
    137   matrix[1] = cell_size[1];
    138   matrix[2] = cell_size[3];
    139   matrix[3] = cell_size[1];
    140   matrix[4] = cell_size[2];
    141   matrix[5] = cell_size[4];
    142   matrix[6] = cell_size[3];
    143   matrix[7] = cell_size[4];
    144   matrix[8] = cell_size[5];
    145   // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    146   for (N[0]=-1;N[0]<=1;N[0]++)
    147     for (N[1]=-1;N[1]<=1;N[1]++)
    148       for (N[2]=-1;N[2]<=1;N[2]++) {
    149         // create the translation vector
    150         TranslationVector.Zero();
    151         for (int i=NDIM;i--;)
    152           TranslationVector.x[i] = (double)N[i];
    153         TranslationVector.MatrixMultiplication(matrix);
    154         // add onto the original vector to compare with
    155         Shiftedy = y + TranslationVector;
    156         // get distance and compare with minimum so far
    157         tmp = DistanceSquared(Shiftedy);
    158         if (tmp < res) res = tmp;
    159       }
    160   return (res);
     101  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     102  return rep->PeriodicDistanceSquared(y,cell_size);
    161103};
    162104
    … …  
    167109void Vector::KeepPeriodic(const double * const matrix)
    168110{
    169 //  int N[NDIM];
    170 //  bool flag = false;
    171   //vector Shifted, TranslationVector;
    172   Vector TestVector;
    173 //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
    174 //  Log() << Verbose(2) << "Vector is: ";
    175 //  Output(out);
    176 //  Log() << Verbose(0) << endl;
    177   TestVector = (*this);
    178   TestVector.InverseMatrixMultiplication(matrix);
    179   for(int i=NDIM;i--;) { // correct periodically
    180     if (TestVector.x[i] < 0) {  // get every coefficient into the interval [0,1)
    181       TestVector.x[i] += ceil(TestVector.x[i]);
    182     } else {
    183       TestVector.x[i] -= floor(TestVector.x[i]);
    184     }
    185   }
    186   TestVector.MatrixMultiplication(matrix);
    187   (*this) = TestVector;
    188 //  Log() << Verbose(2) << "New corrected vector is: ";
    189 //  Output(out);
    190 //  Log() << Verbose(0) << endl;
    191 //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
     111  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     112  rep->KeepPeriodic(matrix);
    192113};
    193114
    … …  
    198119double Vector::ScalarProduct(const Vector &y) const
    199120{
    200   double res = 0.;
    201   for (int i=NDIM;i--;)
    202     res += x[i]*y[i];
    203   return (res);
     121  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     122  return rep->ScalarProduct(y);
    204123};
    205124
    … …  
    213132void Vector::VectorProduct(const Vector &y)
    214133{
    215   Vector tmp;
    216   tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
    217   tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
    218   tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
    219   (*this) = tmp;
     134  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     135  rep->VectorProduct(y);
    220136};
    221137
    … …  
    227143void Vector::ProjectOntoPlane(const Vector &y)
    228144{
    229   Vector tmp = y;
    230   tmp.Normalize();
    231   tmp *= ScalarProduct(tmp);
    232   (*this) -= tmp;
     145  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     146  rep->ProjectOntoPlane(y);
    233147};
    234148
    … …  
    241155double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
    242156{
    243   // first create part that is orthonormal to PlaneNormal with withdraw
    244   Vector temp = (*this) - PlaneOffset;
    245   temp.MakeNormalTo(PlaneNormal);
    246   temp *= -1.;
    247   // then add connecting vector from plane to point
    248   temp += (*this);
    249   temp -= PlaneOffset;
    250   double sign = temp.ScalarProduct(PlaneNormal);
    251   if (fabs(sign) > MYEPSILON)
    252     sign /= fabs(sign);
    253   else
    254     sign = 0.;
    255 
    256   return (temp.Norm()*sign);
     157  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     158  return rep->DistanceToPlane(PlaneNormal,PlaneOffset);
    257159};
    258160
    … …  
    262164void Vector::ProjectIt(const Vector &y)
    263165{
    264   Vector helper = y;
    265   helper.Scale(-(ScalarProduct(y)));
    266   AddVector(helper);
     166  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     167  rep->ProjectIt(y);
    267168};
    268169
    … …  
    273174Vector Vector::Projection(const Vector &y) const
    274175{
    275   Vector helper = y;
    276   helper.Scale((ScalarProduct(y)/y.NormSquared()));
    277 
    278   return helper;
     176  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     177  return rep->Projection(y);
    279178};
    280179
    … …  
    299198void Vector::Normalize()
    300199{
    301   double res = 0.;
    302   for (int i=NDIM;i--;)
    303     res += this->x[i]*this->x[i];
    304   if (fabs(res) > MYEPSILON)
    305     res = 1./sqrt(res);
    306   Scale(&res);
     200  double factor = Norm();
     201  (*this) *= 1/factor;
    307202};
    308203
    … …  
    311206void Vector::Zero()
    312207{
    313   for (int i=NDIM;i--;)
    314     this->x[i] = 0.;
     208  rep.reset(new SingleVector());
    315209};
    316210
    … …  
    319213void Vector::One(const double one)
    320214{
    321   for (int i=NDIM;i--;)
    322     this->x[i] = one;
    323 };
    324 
    325 /** Initialises all components of this vector.
    326  */
    327 void Vector::Init(const double x1, const double x2, const double x3)
    328 {
    329   x[0] = x1;
    330   x[1] = x2;
    331   x[2] = x3;
     215  rep.reset(new SingleVector(one,one,one));
    332216};
    333217
    … …  
    337221bool Vector::IsZero() const
    338222{
    339   return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
     223  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     224  return rep->IsZero();
    340225};
    341226
    … …  
    345230bool Vector::IsOne() const
    346231{
    347   return (fabs(Norm() - 1.) < MYEPSILON);
     232  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     233  return rep->IsOne();
    348234};
    349235
    … …  
    353239bool Vector::IsNormalTo(const Vector &normal) const
    354240{
    355   if (ScalarProduct(normal) < MYEPSILON)
    356     return true;
    357   else
    358     return false;
     241  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     242  return rep->IsNormalTo(normal);
    359243};
    360244
    … …  
    364248bool Vector::IsEqualTo(const Vector &a) const
    365249{
    366   bool status = true;
    367   for (int i=0;i<NDIM;i++) {
    368     if (fabs(x[i] - a[i]) > MYEPSILON)
    369       status = false;
    370   }
    371   return status;
     250  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     251  return rep->IsEqualTo(a);
    372252};
    373253
    … …  
    378258double Vector::Angle(const Vector &y) const
    379259{
    380   double norm1 = Norm(), norm2 = y.Norm();
    381   double angle = -1;
    382   if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
    383     angle = this->ScalarProduct(y)/norm1/norm2;
    384   // -1-MYEPSILON occured due to numerical imprecision, catch ...
    385   //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
    386   if (angle < -1)
    387     angle = -1;
    388   if (angle > 1)
    389     angle = 1;
    390   return acos(angle);
     260  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     261  return rep->Angle(y);
    391262};
    392263
    393264
    394265double& Vector::operator[](size_t i){
    395   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
    396   return x[i];
     266  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     267  return (*rep)[i];
    397268}
    398269
    399270const double& Vector::operator[](size_t i) const{
    400   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
    401   return x[i];
     271  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     272  return (*rep)[i];
    402273}
    403274
    … …  
    411282
    412283double* Vector::get(){
    413   return x;
     284  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     285  return rep->get();
    414286}
    415287
    … …  
    421293bool Vector::operator==(const Vector& b) const
    422294{
    423   bool status = true;
    424   for (int i=0;i<NDIM;i++)
    425     status = status && (fabs((*this)[i] - b[i]) < MYEPSILON);
    426   return status;
     295  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
     296  return IsEqualTo(b);
    427297};
    428298
    … …  
    467337Vector const Vector::operator+(const Vector& b) const
    468338{
     339  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
    469340  Vector x = *this;
    470341  x.AddVector(b);
    … …  
    479350Vector const Vector::operator-(const Vector& b) const
    480351{
     352  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
    481353  Vector x = *this;
    482354  x.SubtractVector(b);
    … …  
    520392};
    521393
    522 /** Scales each atom coordinate by an individual \a factor.
    523  * \param *factor pointer to scaling factor
    524  */
    525 void Vector::Scale(const double ** const factor)
    526 {
    527   for (int i=NDIM;i--;)
    528     x[i] *= (*factor)[i];
    529 };
    530 
    531 void Vector::Scale(const double * const factor)
    532 {
    533   for (int i=NDIM;i--;)
    534     x[i] *= *factor;
    535 };
     394
     395void Vector::ScaleAll(const double *factor)
     396{
     397  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     398  rep->ScaleAll(factor);
     399};
     400
     401
    536402
    537403void Vector::Scale(const double factor)
    538404{
    539   for (int i=NDIM;i--;)
    540     x[i] *= factor;
    541 };
    542 
    543 /** Translate atom by given vector.
    544  * \param trans[] translation vector.
    545  */
    546 void Vector::Translate(const Vector &trans)
    547 {
    548   for (int i=NDIM;i--;)
    549     x[i] += trans[i];
     405  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     406  rep->Scale(factor);
    550407};
    551408
    … …  
    556413void Vector::WrapPeriodically(const double * const M, const double * const Minv)
    557414{
    558   MatrixMultiplication(Minv);
    559   // truncate to [0,1] for each axis
    560   for (int i=0;i<NDIM;i++) {
    561     x[i] += 0.5;  // set to center of box
    562     while (x[i] >= 1.)
    563       x[i] -= 1.;
    564     while (x[i] < 0.)
    565       x[i] += 1.;
    566   }
    567   MatrixMultiplication(M);
     415  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     416  rep->WrapPeriodically(M,Minv);
    568417};
    569418
    … …  
    573422void Vector::MatrixMultiplication(const double * const M)
    574423{
    575   Vector C;
    576   // do the matrix multiplication
    577   C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    578   C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
    579   C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    580   // transfer the result into this
    581   for (int i=NDIM;i--;)
    582     x[i] = C.x[i];
     424  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     425  rep->MatrixMultiplication(M);
    583426};
    584427
    … …  
    588431bool Vector::InverseMatrixMultiplication(const double * const A)
    589432{
    590   Vector C;
    591   double B[NDIM*NDIM];
    592   double detA = RDET3(A);
    593   double detAReci;
    594 
    595   // calculate the inverse B
    596   if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    597     detAReci = 1./detA;
    598     B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    599     B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    600     B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    601     B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    602     B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    603     B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    604     B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    605     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    606     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
    607 
    608     // do the matrix multiplication
    609     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
    610     C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
    611     C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    612     // transfer the result into this
    613     for (int i=NDIM;i--;)
    614       x[i] = C.x[i];
    615     return true;
    616   } else {
    617     return false;
    618   }
     433  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     434  return rep->InverseMatrixMultiplication(A);
    619435};
    620436
    … …  
    639455void Vector::Mirror(const Vector &n)
    640456{
    641   double projection;
    642   projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
    643   // withdraw projected vector twice from original one
    644   for (int i=NDIM;i--;)
    645     x[i] -= 2.*projection*n[i];
     457  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     458  rep->Mirror(n);
    646459};
    647460
    … …  
    655468bool Vector::MakeNormalTo(const Vector &y1)
    656469{
    657   bool result = false;
    658   double factor = y1.ScalarProduct(*this)/y1.NormSquared();
    659   Vector x1 = factor * y1 ;
    660   SubtractVector(x1);
    661   for (int i=NDIM;i--;)
    662     result = result || (fabs(x[i]) > MYEPSILON);
    663 
    664   return result;
     470  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     471  return rep->MakeNormalTo(y1);
    665472};
    666473
    … …  
    673480bool Vector::GetOneNormalVector(const Vector &GivenVector)
    674481{
    675   int Components[NDIM]; // contains indices of non-zero components
    676   int Last = 0;   // count the number of non-zero entries in vector
    677   int j;  // loop variables
    678   double norm;
    679 
    680   for (j=NDIM;j--;)
    681     Components[j] = -1;
    682   // find two components != 0
    683   for (j=0;j<NDIM;j++)
    684     if (fabs(GivenVector[j]) > MYEPSILON)
    685       Components[Last++] = j;
    686 
    687   switch(Last) {
    688     case 3:  // threecomponent system
    689     case 2:  // two component system
    690       norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
    691       x[Components[2]] = 0.;
    692       // in skp both remaining parts shall become zero but with opposite sign and third is zero
    693       x[Components[1]] = -1./GivenVector[Components[1]] / norm;
    694       x[Components[0]] = 1./GivenVector[Components[0]] / norm;
    695       return true;
    696       break;
    697     case 1: // one component system
    698       // set sole non-zero component to 0, and one of the other zero component pendants to 1
    699       x[(Components[0]+2)%NDIM] = 0.;
    700       x[(Components[0]+1)%NDIM] = 1.;
    701       x[Components[0]] = 0.;
    702       return true;
    703       break;
    704     default:
    705       return false;
    706   }
     482  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     483  return rep->GetOneNormalVector(GivenVector);
    707484};
    708485
    … …  
    712489void Vector::AddVector(const Vector &y)
    713490{
    714   for (int i=NDIM;i--;)
    715     this->x[i] += y[i];
     491  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     492  rep->AddVector(y);
    716493}
    717494
    … …  
    721498void Vector::SubtractVector(const Vector &y)
    722499{
    723   for (int i=NDIM;i--;)
    724     this->x[i] -= y[i];
    725 }
    726 
    727 /** Copy vector \a y componentwise.
    728  * \param y vector
    729  */
    730 void Vector::CopyVector(const Vector &y)
    731 {
    732   // check for self assignment
    733   if(&y!=this) {
    734     for (int i=NDIM;i--;)
    735       this->x[i] = y.x[i];
    736   }
    737 }
    738 
    739 // this function is completely unused so it is deactivated until new uses arrive and a new
    740 // place can be found
    741 #if 0
    742 /** Solves a vectorial system consisting of two orthogonal statements and a norm statement.
    743  * This is linear system of equations to be solved, however of the three given (skp of this vector\
    744  * with either of the three hast to be zero) only two are linear independent. The third equation
    745  * is that the vector should be of magnitude 1 (orthonormal). This all leads to a case-based solution
    746  * where very often it has to be checked whether a certain value is zero or not and thus forked into
    747  * another case.
    748  * \param *x1 first vector
    749  * \param *x2 second vector
    750  * \param *y third vector
    751  * \param alpha first angle
    752  * \param beta second angle
    753  * \param c norm of final vector
    754  * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
    755  * \bug this is not yet working properly
    756  */
    757 bool Vector::SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c)
    758 {
    759   double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    760   double ang; // angle on testing
    761   double sign[3];
    762   int i,j,k;
    763   A = cos(alpha) * x1->Norm() * c;
    764   B1 = cos(beta + M_PI/2.) * y->Norm() * c;
    765   B2 = cos(beta) * x2->Norm() * c;
    766   C = c * c;
    767   Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
    768   int flag = 0;
    769   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    770     if (fabs(x1->x[1]) > MYEPSILON) {
    771       flag = 1;
    772     } else if (fabs(x1->x[2]) > MYEPSILON) {
    773        flag = 2;
    774     } else {
    775       return false;
    776     }
    777   }
    778   switch (flag) {
    779     default:
    780     case 0:
    781       break;
    782     case 2:
    783       flip(x1->x[0],x1->x[1]);
    784       flip(x2->x[0],x2->x[1]);
    785       flip(y->x[0],y->x[1]);
    786       //flip(x[0],x[1]);
    787       flip(x1->x[1],x1->x[2]);
    788       flip(x2->x[1],x2->x[2]);
    789       flip(y->x[1],y->x[2]);
    790       //flip(x[1],x[2]);
    791     case 1:
    792       flip(x1->x[0],x1->x[1]);
    793       flip(x2->x[0],x2->x[1]);
    794       flip(y->x[0],y->x[1]);
    795       //flip(x[0],x[1]);
    796       flip(x1->x[1],x1->x[2]);
    797       flip(x2->x[1],x2->x[2]);
    798       flip(y->x[1],y->x[2]);
    799       //flip(x[1],x[2]);
    800       break;
    801   }
    802   // now comes the case system
    803   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
    804   D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    805   D3 = y->x[0]/x1->x[0]*A-B1;
    806   Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
    807   if (fabs(D1) < MYEPSILON) {
    808     Log() << Verbose(2) << "D1 == 0!\n";
    809     if (fabs(D2) > MYEPSILON) {
    810       Log() << Verbose(3) << "D2 != 0!\n";
    811       x[2] = -D3/D2;
    812       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    813       E2 = -x1->x[1]/x1->x[0];
    814       Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    815       F1 = E1*E1 + 1.;
    816       F2 = -E1*E2;
    817       F3 = E1*E1 + D3*D3/(D2*D2) - C;
    818       Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    819       if (fabs(F1) < MYEPSILON) {
    820         Log() << Verbose(4) << "F1 == 0!\n";
    821         Log() << Verbose(4) << "Gleichungssystem linear\n";
    822         x[1] = F3/(2.*F2);
    823       } else {
    824         p = F2/F1;
    825         q = p*p - F3/F1;
    826         Log() << Verbose(4) << "p " << p << "\tq " << q << endl;
    827         if (q < 0) {
    828           Log() << Verbose(4) << "q < 0" << endl;
    829           return false;
    830         }
    831         x[1] = p + sqrt(q);
    832       }
    833       x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    834     } else {
    835       Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n";
    836       return false;
    837     }
    838   } else {
    839     E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    840     E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    841     Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    842     F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    843     F2 = -(E1*E2 + D2*D3/(D1*D1));
    844     F3 = E1*E1 + D3*D3/(D1*D1) - C;
    845     Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    846     if (fabs(F1) < MYEPSILON) {
    847       Log() << Verbose(3) << "F1 == 0!\n";
    848       Log() << Verbose(3) << "Gleichungssystem linear\n";
    849       x[2] = F3/(2.*F2);
    850     } else {
    851       p = F2/F1;
    852       q = p*p - F3/F1;
    853       Log() << Verbose(3) << "p " << p << "\tq " << q << endl;
    854       if (q < 0) {
    855         Log() << Verbose(3) << "q < 0" << endl;
    856         return false;
    857       }
    858       x[2] = p + sqrt(q);
    859     }
    860     x[1] = (-D2 * x[2] - D3)/D1;
    861     x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    862   }
    863   switch (flag) { // back-flipping
    864     default:
    865     case 0:
    866       break;
    867     case 2:
    868       flip(x1->x[0],x1->x[1]);
    869       flip(x2->x[0],x2->x[1]);
    870       flip(y->x[0],y->x[1]);
    871       flip(x[0],x[1]);
    872       flip(x1->x[1],x1->x[2]);
    873       flip(x2->x[1],x2->x[2]);
    874       flip(y->x[1],y->x[2]);
    875       flip(x[1],x[2]);
    876     case 1:
    877       flip(x1->x[0],x1->x[1]);
    878       flip(x2->x[0],x2->x[1]);
    879       flip(y->x[0],y->x[1]);
    880       //flip(x[0],x[1]);
    881       flip(x1->x[1],x1->x[2]);
    882       flip(x2->x[1],x2->x[2]);
    883       flip(y->x[1],y->x[2]);
    884       flip(x[1],x[2]);
    885       break;
    886   }
    887   // one z component is only determined by its radius (without sign)
    888   // thus check eight possible sign flips and determine by checking angle with second vector
    889   for (i=0;i<8;i++) {
    890     // set sign vector accordingly
    891     for (j=2;j>=0;j--) {
    892       k = (i & pot(2,j)) << j;
    893       Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
    894       sign[j] = (k == 0) ? 1. : -1.;
    895     }
    896     Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
    897     // apply sign matrix
    898     for (j=NDIM;j--;)
    899       x[j] *= sign[j];
    900     // calculate angle and check
    901     ang = x2->Angle (this);
    902     Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
    903     if (fabs(ang - cos(beta)) < MYEPSILON) {
    904       break;
    905     }
    906     // unapply sign matrix (is its own inverse)
    907     for (j=NDIM;j--;)
    908       x[j] *= sign[j];
    909   }
    910   return true;
    911 };
    912 
    913 #endif
     500  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     501  rep->SubtractVector(y);
     502}
    914503
    915504/**
    … …  
    922511bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
    923512{
    924   Vector a = (*this) - offset;
    925   a.InverseMatrixMultiplication(parallelepiped);
    926   bool isInside = true;
    927 
    928   for (int i=NDIM;i--;)
    929     isInside = isInside && ((a.x[i] <= 1) && (a.x[i] >= 0));
    930 
    931   return isInside;
    932 }
     513  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
     514  return rep->IsInParallelepiped(offset, parallelepiped);
     515}
     516
     517bool Vector::isBaseClass() const{
     518  return true;
     519}
     520
     521Vector* Vector::clone() const{
     522  ASSERT(false, "Cannot clone a base Vector object");
     523  return 0;
     524}

  • src/vector.hpp

    r273382 r1bd79e  
    1515#include <gsl/gsl_multimin.h>
    1616
     17#include <memory>
     18
    1719#include "defs.hpp"
    1820
    … …  
    2628  // this struct is used to indicate calls to the Baseconstructor from inside vectors.
    2729  struct Baseconstructor{};
    28   public:
     30public:
    2931
    3032  Vector();
    … …  
    3335  virtual ~Vector();
    3436
    35   virtual double Distance(const Vector &y) const;
     37  // Method implemented by forwarding to the Representation
     38
    3639  virtual double DistanceSquared(const Vector &y) const;
    3740  virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
    … …  
    3942  virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
    4043  virtual double ScalarProduct(const Vector &y) const;
    41   virtual double Norm() const;
    42   virtual double NormSquared() const;
    4344  virtual double Angle(const Vector &y) const;
    4445  virtual bool IsZero() const;
    … …  
    4950  virtual void AddVector(const Vector &y);
    5051  virtual void SubtractVector(const Vector &y);
    51   virtual void CopyVector(const Vector &y);
    5252  virtual void VectorProduct(const Vector &y);
    5353  virtual void ProjectOntoPlane(const Vector &y);
    5454  virtual void ProjectIt(const Vector &y);
    5555  virtual Vector Projection(const Vector &y) const;
    56   virtual void Zero();
    57   virtual void One(const double one);
    58   virtual void Init(const double x1, const double x2, const double x3);
    59   virtual void Normalize();
    60   virtual void Translate(const Vector &x);
    6156  virtual void Mirror(const Vector &x);
    62   virtual void Scale(const double ** const factor);
    63   virtual void Scale(const double * const factor);
     57  virtual void ScaleAll(const double *factor);
    6458  virtual void Scale(const double factor);
    6559  virtual void MatrixMultiplication(const double * const M);
    6660  virtual bool InverseMatrixMultiplication(const double * const M);
    6761  virtual void KeepPeriodic(const double * const matrix);
    68   virtual void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
    6962  virtual bool GetOneNormalVector(const Vector &x1);
    7063  virtual bool MakeNormalTo(const Vector &y1);
    71   //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
    7264  virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
    7365  virtual void WrapPeriodically(const double * const M, const double * const Minv);
    … …  
    7668  virtual double& operator[](size_t i);
    7769  virtual const double& operator[](size_t i) const;
    78   virtual double& at(size_t i);
    79   virtual const double& at(size_t i) const;
     70  double& at(size_t i);
     71  const double& at(size_t i) const;
    8072
    8173  // Assignment operator
    … …  
    8476  // Access to internal structure
    8577  virtual double* get();
     78
     79  // Methods that are derived directly from other methods
     80  double Distance(const Vector &y) const;
     81  double Norm() const;
     82  double NormSquared() const;
     83  void Normalize();
     84  void Zero();
     85  void One(const double one);
     86  void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
    8687
    8788  // operators for mathematical operations
    … …  
    9293  Vector const operator-(const Vector& b) const;
    9394
     95protected:
     96  typedef std::auto_ptr<Vector> rep_ptr;
     97  Vector(Baseconstructor);
     98  Vector(Baseconstructor,const Vector*);
     99  static Vector VecFromRep(const Vector*);
     100
    94101private:
    95   double x[NDIM];
     102  // method used for protection, i.e. to avoid infinite recursion
     103  // when our internal rep becomes messed up
     104  virtual bool isBaseClass() const;
     105  virtual Vector* clone() const;
     106  // this is used to represent the vector internally
     107  rep_ptr rep;
    96108
    97109};

  • src/vector_ops.cpp

    r273382 r1bd79e  
    120120  Vector res;
    121121  // normalise this vector with respect to axis
    122   a.CopyVector(vec);
     122  a = vec;
    123123  a.ProjectOntoPlane(axis);
    124124  // construct normal vector
Note: See TracChangeset for help on using the changeset viewer.