Changes in src/Box.cpp [57f243:0ae36b]

File:

• src/Box.cpp (modified) (8 diffs)

src/Box.cpp

@@ -11 +11 @@
 
 #include <cmath>
-#include <iostream>
-#include <cstdlib>
 
-#include "LinearAlgebra/Matrix.hpp"
-#include "LinearAlgebra/Vector.hpp"
-#include "LinearAlgebra/Plane.hpp"
+#include "Matrix.hpp"
+#include "vector.hpp"
 
 #include "Helpers/Assert.hpp"
-
-using namespace std;
 
 Box::Box()
@@ -28 +23 @@
   Minv = new Matrix();
   Minv->one();
-  conditions.resize(3);
-  conditions[0] = conditions[1] = conditions[2] = Wrap;
 }
 
@@ -35 +28 @@
   M=new Matrix(*src.M);
   Minv = new Matrix(*src.Minv);
-  conditions = src.conditions;
 }
 
@@ -68 +60 @@
   Vector helper = translateOut(point);
   for(int i=NDIM;i--;){
-
-    switch(conditions[i]){
-    case Wrap:
-      helper.at(i)=fmod(helper.at(i),1);
-      helper.at(i)+=(helper.at(i)>=0)?0:1;
-      break;
-    case Bounce:
-      {
-        // there probably is a better way to handle this...
-        // all the fabs and fmod modf probably makes it very slow
-        double intpart,fracpart;
-        fracpart = modf(fabs(helper.at(i)),&intpart);
-        helper.at(i) = fabs(fracpart-fmod(intpart,2));
-      }
-      break;
-    case Ignore:
-      break;
-    default:
-      ASSERT(0,"No default case for this");
-    }
-
+    double intpart,fracpart;
+    fracpart = modf(helper.at(i),&intpart);
+    if(fracpart<0.)
+      fracpart+=1.;
+    helper.at(i)=fracpart;
   }
   return translateIn(helper);
@@ -96 +72 @@
 {
   bool result = true;
-  Vector tester = translateOut(point);
+  Vector tester = (*Minv) * point;
 
   for(int i=0;i<NDIM;i++)
-    result = result &&
-              ((conditions[i] == Ignore) ||
-               ((tester[i] >= -MYEPSILON) &&
-               ((tester[i] - 1.) < MYEPSILON)));
+    result = result && ((tester[i] >= -MYEPSILON) && ((tester[i] - 1.) < MYEPSILON));
 
   return result;
@@ -109 +82 @@
 
 VECTORSET(std::list) Box::explode(const Vector &point,int n) const{
-  ASSERT(isInside(point),"Exploded point not inside Box");
   VECTORSET(std::list) res;
 
-  Vector translater = translateOut(point);
-  Vector mask; // contains the ignored coordinates
+  // translate the Vector into each of the 27 neighbourhoods
 
-  // count the number of coordinates we need to do
-  int dims = 0; // number of dimensions that are not ignored
-  vector<int> coords;
-  vector<int> index;
-  for(int i=0;i<NDIM;++i){
-    if(conditions[i]==Ignore){
-      mask[i]=translater[i];
-      continue;
-    }
-    coords.push_back(i);
-    index.push_back(-n);
-    dims++;
-  } // there are max vectors in total we need to create
-
-  if(!dims){
-    // all boundaries are ignored
-    res.push_back(point);
-    return res;
+  // first create all translation Vectors
+  // there are (n*2+1)^3 such vectors
+  int max_dim = (n*2+1);
+  int max_dim2 = max_dim*max_dim;
+  int max = max_dim2*max_dim;
+  // only one loop to avoid unneccessary jumps
+  for(int i = 0;i<max;++i){
+    // get all coordinates for this iteration
+    int n1 = (i%max_dim)-n;
+    int n2 = ((i/max_dim)%max_dim)-n;
+    int n3 = ((i/max_dim2))-n;
+    Vector translation = translateIn(Vector(n1,n2,n3));
+    res.push_back(translation);
   }
-
-  bool done = false;
-  while(!done){
-    // create this vector
-    Vector helper;
-    for(int i=0;i<dims;++i){
-      switch(conditions[coords[i]]){
-        case Wrap:
-          helper[coords[i]] = index[i]+translater[coords[i]];
-          break;
-        case Bounce:
-          {
-            // Bouncing the coordinate x produces the series:
-            // 0 -> x
-            // 1 -> 2-x
-            // 2 -> 2+x
-            // 3 -> 4-x
-            // 4 -> 4+x
-            // the first number is the next bigger even number (n+n%2)
-            // the next number is the value with alternating sign (x-2*(n%2)*x)
-            // the negative numbers produce the same sequence reversed and shifted
-            int n = abs(index[i]) + ((index[i]<0)?-1:0);
-            int sign = (index[i]<0)?-1:+1;
-            int even = n%2;
-            helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
-            helper[coords[i]]*=sign;
-          }
-          break;
-        case Ignore:
-          ASSERT(0,"Ignored coordinate handled in generation loop");
-        default:
-          ASSERT(0,"No default case for this switch-case");
-      }
-
-    }
-    // add back all ignored coordinates (not handled in above loop)
-    helper+=mask;
-    res.push_back(translateIn(helper));
-    // set the new indexes
-    int pos=0;
-    ++index[pos];
-    while(index[pos]>n){
-      index[pos++]=-n;
-      if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
-        done = true;
-        break;
-      }
-      ++index[pos];
-    }
-  }
+  // translate all the translation vector by the offset defined by point
+  res.translate(point);
   return res;
 }
 
 VECTORSET(std::list) Box::explode(const Vector &point) const{
-  ASSERT(isInside(point),"Exploded point not inside Box");
-  return explode(point,1);
+  VECTORSET(std::list) res;
+
+  // translate the Vector into each of the 27 neighbourhoods
+
+  // first create all 27 translation Vectors
+  // these loops depend on fixed parameters and can easily be expanded
+  // by the compiler to allow code without jumps
+  for(int n1 = -1;n1<=1;++n1){
+    for(int n2 = -1;n2<=1;++n2){
+      for(int n3 = -1;n3<=1;++n3){
+        // get all coordinates for this iteration
+        Vector translation = translateIn(Vector(n1,n2,n3));
+        res.push_back(translation);
+      }
+    }
+  }
+  // translate all the translation vector by the offset defined by point
+  res.translate(point);
+  return res;
 }
 
 double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
-  Vector helper1 = WrapPeriodically(point1);
-  Vector helper2 = WrapPeriodically(point2);
-  VECTORSET(std::list) expansion = explode(helper1);
-  double res = expansion.minDistSquared(helper2);
+  VECTORSET(std::list) expansion = explode(point1);
+  double res = expansion.minDistSquared(point2);
   return res;
 }
@@ -206 +139 @@
 }
 
-const Box::Conditions_t Box::getConditions(){
-  return conditions;
-}
-
-void Box::setCondition(int i,Box::BoundaryCondition_t condition){
-  conditions[i]=condition;
-}
-
-const vector<pair<Plane,Plane> >  Box::getBoundingPlanes(){
-  vector<pair<Plane,Plane> > res;
-  for(int i=0;i<NDIM;++i){
-    Vector base1,base2,base3;
-    base2[(i+1)%NDIM] = 1.;
-    base3[(i+2)%NDIM] = 1.;
-    Plane p1(translateIn(base1),
-             translateIn(base2),
-             translateIn(base3));
-    Vector offset;
-    offset[i]=1;
-    Plane p2(translateIn(base1+offset),
-             translateIn(base2+offset),
-             translateIn(base3+offset));
-    res.push_back(make_pair(p1,p2));
-  }
-  return res;
-}
-
-void Box::setCuboid(const Vector &endpoint){
-  ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
-  M->one();
-  M->diagonal()=endpoint;
-  Vector &dinv = Minv->diagonal();
-  for(int i=NDIM;i--;)
-    dinv[i]=1/endpoint[i];
-}
-
 Box &Box::operator=(const Box &src){
   if(&src!=this){
@@ -248 +145 @@
     M    = new Matrix(*src.M);
     Minv = new Matrix(*src.Minv);
-    conditions = src.conditions;
   }
   return *this;