Changeset 753f02 for src/vector.cpp

Timestamp:

Apr 28, 2010, 2:52:56 PM (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:

8cbb97

Parents:

1bd79e

Message:

Removed Algebraic Hierachy from vector.

File:

• src/vector.cpp (modified) (35 diffs)

src/vector.cpp

@@ -6 +6 @@
 
 
-#include "SingleVector.hpp"
+#include "vector.hpp"
 #include "Helpers/Assert.hpp"
+#include "Helpers/fast_functions.hpp"
 
 #include <iostream>
@@ -18 +19 @@
 /** Constructor of class vector.
  */
-Vector::Vector() :
-  rep(new SingleVector())
-{};
-
-Vector::Vector(Baseconstructor)  // used by derived objects to construct their bases
-{}
-
-Vector::Vector(Baseconstructor,const Vector* v) :
-  rep(v->clone())
-{}
-
-Vector Vector::VecFromRep(const Vector* v){
-  return Vector(Baseconstructor(),v);
-}
-
-/** Constructor of class vector.
- */
-Vector::Vector(const double x1, const double x2, const double x3) :
-  rep(new SingleVector(x1,x2,x3))
-{};
+Vector::Vector()
+{
+  x[0] = x[1] = x[2] = 0.;
+};
 
 /**
  * Copy constructor
  */
-Vector::Vector(const Vector& src) :
-  rep(src.rep->clone())
-{}
+
+Vector::Vector(const Vector& src)
+{
+  x[0] = src[0];
+  x[1] = src[1];
+  x[2] = src[2];
+}
+
+/** Constructor of class vector.
+ */
+Vector::Vector(const double x1, const double x2, const double x3)
+{
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
 
 /**
@@ -50 +48 @@
  */
 Vector& Vector::operator=(const Vector& src){
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   // check for self assignment
   if(&src!=this){
-    rep.reset(src.rep->clone());
+    x[0] = src[0];
+    x[1] = src[1];
+    x[2] = src[2];
   }
   return *this;
@@ -68 +67 @@
 double Vector::DistanceSquared(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->DistanceSquared(y);
+  double res = 0.;
+  for (int i=NDIM;i--;)
+    res += (x[i]-y[i])*(x[i]-y[i]);
+  return (res);
 };
 
@@ -88 +89 @@
 double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->PeriodicDistance(y,cell_size);
+  double res = Distance(y), tmp, matrix[NDIM*NDIM];
+    Vector Shiftedy, TranslationVector;
+    int N[NDIM];
+    matrix[0] = cell_size[0];
+    matrix[1] = cell_size[1];
+    matrix[2] = cell_size[3];
+    matrix[3] = cell_size[1];
+    matrix[4] = cell_size[2];
+    matrix[5] = cell_size[4];
+    matrix[6] = cell_size[3];
+    matrix[7] = cell_size[4];
+    matrix[8] = cell_size[5];
+    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
+    for (N[0]=-1;N[0]<=1;N[0]++)
+      for (N[1]=-1;N[1]<=1;N[1]++)
+        for (N[2]=-1;N[2]<=1;N[2]++) {
+          // create the translation vector
+          TranslationVector.Zero();
+          for (int i=NDIM;i--;)
+            TranslationVector[i] = (double)N[i];
+          TranslationVector.MatrixMultiplication(matrix);
+          // add onto the original vector to compare with
+          Shiftedy = y + TranslationVector;
+          // get distance and compare with minimum so far
+          tmp = Distance(Shiftedy);
+          if (tmp < res) res = tmp;
+        }
+    return (res);
 };
 
@@ -99 +126 @@
 double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->PeriodicDistanceSquared(y,cell_size);
+  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
+    Vector Shiftedy, TranslationVector;
+    int N[NDIM];
+    matrix[0] = cell_size[0];
+    matrix[1] = cell_size[1];
+    matrix[2] = cell_size[3];
+    matrix[3] = cell_size[1];
+    matrix[4] = cell_size[2];
+    matrix[5] = cell_size[4];
+    matrix[6] = cell_size[3];
+    matrix[7] = cell_size[4];
+    matrix[8] = cell_size[5];
+    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
+    for (N[0]=-1;N[0]<=1;N[0]++)
+      for (N[1]=-1;N[1]<=1;N[1]++)
+        for (N[2]=-1;N[2]<=1;N[2]++) {
+          // create the translation vector
+          TranslationVector.Zero();
+          for (int i=NDIM;i--;)
+            TranslationVector[i] = (double)N[i];
+          TranslationVector.MatrixMultiplication(matrix);
+          // add onto the original vector to compare with
+          Shiftedy = y + TranslationVector;
+          // get distance and compare with minimum so far
+          tmp = DistanceSquared(Shiftedy);
+          if (tmp < res) res = tmp;
+        }
+    return (res);
 };
 
@@ -109 +162 @@
 void Vector::KeepPeriodic(const double * const matrix)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->KeepPeriodic(matrix);
+  //  int N[NDIM];
+  //  bool flag = false;
+    //vector Shifted, TranslationVector;
+  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
+  //  Log() << Verbose(2) << "Vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+    InverseMatrixMultiplication(matrix);
+    for(int i=NDIM;i--;) { // correct periodically
+      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
+        at(i) += ceil(at(i));
+      } else {
+        at(i) -= floor(at(i));
+      }
+    }
+    MatrixMultiplication(matrix);
+  //  Log() << Verbose(2) << "New corrected vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
 };
 
@@ -119 +190 @@
 double Vector::ScalarProduct(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->ScalarProduct(y);
+  double res = 0.;
+  for (int i=NDIM;i--;)
+    res += x[i]*y[i];
+  return (res);
 };
 
@@ -132 +205 @@
 void Vector::VectorProduct(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->VectorProduct(y);
+  Vector tmp;
+  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
+  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
+  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
+  (*this) = tmp;
 };
 
@@ -143 +219 @@
 void Vector::ProjectOntoPlane(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ProjectOntoPlane(y);
+  Vector tmp;
+  tmp = y;
+  tmp.Normalize();
+  tmp.Scale(ScalarProduct(tmp));
+  *this -= tmp;
 };
 
@@ -155 +234 @@
 double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->DistanceToPlane(PlaneNormal,PlaneOffset);
+  // first create part that is orthonormal to PlaneNormal with withdraw
+  Vector temp = (*this) - PlaneOffset;
+  temp.MakeNormalTo(PlaneNormal);
+  temp.Scale(-1.);
+  // then add connecting vector from plane to point
+  temp += (*this)-PlaneOffset;
+  double sign = temp.ScalarProduct(PlaneNormal);
+  if (fabs(sign) > MYEPSILON)
+    sign /= fabs(sign);
+  else
+    sign = 0.;
+
+  return (temp.Norm()*sign);
 };
 
@@ -164 +254 @@
 void Vector::ProjectIt(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ProjectIt(y);
+  (*this) += (-ScalarProduct(y))*y;
 };
 
@@ -174 +263 @@
 Vector Vector::Projection(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->Projection(y);
+  Vector helper = y;
+  helper.Scale((ScalarProduct(y)/y.NormSquared()));
+
+  return helper;
 };
 
@@ -206 +297 @@
 void Vector::Zero()
 {
-  rep.reset(new SingleVector());
+  at(0)=at(1)=at(2)=0;
 };
 
@@ -213 +304 @@
 void Vector::One(const double one)
 {
-  rep.reset(new SingleVector(one,one,one));
+  at(0)=at(1)=at(2)=one;
 };
 
@@ -221 +312 @@
 bool Vector::IsZero() const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsZero();
+  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
 };
 
@@ -230 +320 @@
 bool Vector::IsOne() const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsOne();
+  return (fabs(Norm() - 1.) < MYEPSILON);
 };
 
@@ -239 +328 @@
 bool Vector::IsNormalTo(const Vector &normal) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsNormalTo(normal);
+  if (ScalarProduct(normal) < MYEPSILON)
+    return true;
+  else
+    return false;
 };
 
@@ -248 +339 @@
 bool Vector::IsEqualTo(const Vector &a) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsEqualTo(a);
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };
 
@@ -258 +353 @@
 double Vector::Angle(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->Angle(y);
+  double norm1 = Norm(), norm2 = y.Norm();
+  double angle = -1;
+  if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
+    angle = this->ScalarProduct(y)/norm1/norm2;
+  // -1-MYEPSILON occured due to numerical imprecision, catch ...
+  //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
+  if (angle < -1)
+    angle = -1;
+  if (angle > 1)
+    angle = 1;
+  return acos(angle);
 };
 
 
 double& Vector::operator[](size_t i){
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return (*rep)[i];
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
 }
 
 const double& Vector::operator[](size_t i) const{
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return (*rep)[i];
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
 }
 
@@ -282 +386 @@
 
 double* Vector::get(){
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->get();
+  return x;
 }
 
@@ -293 +396 @@
 bool Vector::operator==(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   return IsEqualTo(b);
 };
@@ -337 +439 @@
 Vector const Vector::operator+(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.AddVector(b);
@@ -350 +451 @@
 Vector const Vector::operator-(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.SubtractVector(b);
@@ -395 +495 @@
 void Vector::ScaleAll(const double *factor)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ScaleAll(factor);
+  for (int i=NDIM;i--;)
+    x[i] *= factor[i];
 };
 
@@ -403 +503 @@
 void Vector::Scale(const double factor)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->Scale(factor);
+  for (int i=NDIM;i--;)
+    x[i] *= factor;
 };
 
@@ -413 +513 @@
 void Vector::WrapPeriodically(const double * const M, const double * const Minv)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->WrapPeriodically(M,Minv);
+  MatrixMultiplication(Minv);
+  // truncate to [0,1] for each axis
+  for (int i=0;i<NDIM;i++) {
+    x[i] += 0.5;  // set to center of box
+    while (x[i] >= 1.)
+      x[i] -= 1.;
+    while (x[i] < 0.)
+      x[i] += 1.;
+  }
+  MatrixMultiplication(M);
 };
 
@@ -422 +530 @@
 void Vector::MatrixMultiplication(const double * const M)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->MatrixMultiplication(M);
+  // do the matrix multiplication
+  at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
+  at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
+  at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
 };
 
@@ -431 +541 @@
 bool Vector::InverseMatrixMultiplication(const double * const A)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->InverseMatrixMultiplication(A);
+  double B[NDIM*NDIM];
+  double detA = RDET3(A);
+  double detAReci;
+
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+
+    // do the matrix multiplication
+    at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
+    at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
+    at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
+
+    return true;
+  } else {
+    return false;
+  }
 };
 
@@ -455 +589 @@
 void Vector::Mirror(const Vector &n)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->Mirror(n);
+  double projection;
+  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
+  // withdraw projected vector twice from original one
+  for (int i=NDIM;i--;)
+    x[i] -= 2.*projection*n[i];
 };
 
@@ -468 +605 @@
 bool Vector::MakeNormalTo(const Vector &y1)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->MakeNormalTo(y1);
+  bool result = false;
+  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
+  Vector x1;
+  x1 = factor * y1;
+  SubtractVector(x1);
+  for (int i=NDIM;i--;)
+    result = result || (fabs(x[i]) > MYEPSILON);
+
+  return result;
 };
 
@@ -480 +624 @@
 bool Vector::GetOneNormalVector(const Vector &GivenVector)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->GetOneNormalVector(GivenVector);
+  int Components[NDIM]; // contains indices of non-zero components
+  int Last = 0;   // count the number of non-zero entries in vector
+  int j;  // loop variables
+  double norm;
+
+  for (j=NDIM;j--;)
+    Components[j] = -1;
+  // find two components != 0
+  for (j=0;j<NDIM;j++)
+    if (fabs(GivenVector[j]) > MYEPSILON)
+      Components[Last++] = j;
+
+  switch(Last) {
+    case 3:  // threecomponent system
+    case 2:  // two component system
+      norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
+      x[Components[2]] = 0.;
+      // in skp both remaining parts shall become zero but with opposite sign and third is zero
+      x[Components[1]] = -1./GivenVector[Components[1]] / norm;
+      x[Components[0]] = 1./GivenVector[Components[0]] / norm;
+      return true;
+      break;
+    case 1: // one component system
+      // set sole non-zero component to 0, and one of the other zero component pendants to 1
+      x[(Components[0]+2)%NDIM] = 0.;
+      x[(Components[0]+1)%NDIM] = 1.;
+      x[Components[0]] = 0.;
+      return true;
+      break;
+    default:
+      return false;
+  }
 };
 
@@ -489 +663 @@
 void Vector::AddVector(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->AddVector(y);
+  for(int i=NDIM;i--;)
+    x[i] += y[i];
 }
 
@@ -498 +672 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->SubtractVector(y);
+  for(int i=NDIM;i--;)
+    x[i] -= y[i];
 }
 
@@ -511 +685 @@
 bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsInParallelepiped(offset, parallelepiped);
-}
-
-bool Vector::isBaseClass() const{
-  return true;
-}
-
-Vector* Vector::clone() const{
-  ASSERT(false, "Cannot clone a base Vector object");
-  return 0;
-}
+  Vector a = (*this)-offset;
+  a.InverseMatrixMultiplication(parallelepiped);
+  bool isInside = true;
+
+  for (int i=NDIM;i--;)
+    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
+
+  return isInside;
+}