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

File:

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

src/Box.cpp

@@ -11 +11 @@
 
 #include <cmath>
-
-#include "Matrix.hpp"
-#include "vector.hpp"
+#include <iostream>
+#include <cstdlib>
+
+#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Plane.hpp"
 
 #include "Helpers/Assert.hpp"
+
+using namespace std;
 
 Box::Box()
@@ -23 +28 @@
   Minv = new Matrix();
   Minv->one();
+  conditions.resize(3);
+  conditions[0] = conditions[1] = conditions[2] = Wrap;
 }
 
@@ -28 +35 @@
   M=new Matrix(*src.M);
   Minv = new Matrix(*src.Minv);
+  conditions = src.conditions;
 }
 
@@ -60 +68 @@
   Vector helper = translateOut(point);
   for(int i=NDIM;i--;){
-    double intpart,fracpart;
-    fracpart = modf(helper.at(i),&intpart);
-    if(fracpart<0.)
-      fracpart+=1.;
-    helper.at(i)=fracpart;
+
+    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");
+    }
+
   }
   return translateIn(helper);
@@ -72 +96 @@
 {
   bool result = true;
-  Vector tester = (*Minv) * point;
+  Vector tester = translateOut(point);
 
   for(int i=0;i<NDIM;i++)
-    result = result && ((tester[i] >= -MYEPSILON) && ((tester[i] - 1.) < MYEPSILON));
+    result = result &&
+              ((conditions[i] == Ignore) ||
+               ((tester[i] >= -MYEPSILON) &&
+               ((tester[i] - 1.) < MYEPSILON)));
 
   return result;
@@ -82 +109 @@
 
 VECTORSET(std::list) Box::explode(const Vector &point,int n) const{
+  ASSERT(isInside(point),"Exploded point not inside Box");
   VECTORSET(std::list) res;
 
-  // translate the Vector into each of the 27 neighbourhoods
-
-  // 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);
-  }
-  // translate all the translation vector by the offset defined by point
-  res.translate(point);
+  Vector translater = translateOut(point);
+  Vector mask; // contains the ignored coordinates
+
+  // 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;
+  }
+
+  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];
+    }
+  }
   return res;
 }
 
 VECTORSET(std::list) Box::explode(const Vector &point) const{
-  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;
+  ASSERT(isInside(point),"Exploded point not inside Box");
+  return explode(point,1);
 }
 
 double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
-  VECTORSET(std::list) expansion = explode(point1);
-  double res = expansion.minDistSquared(point2);
+  Vector helper1 = WrapPeriodically(point1);
+  Vector helper2 = WrapPeriodically(point2);
+  VECTORSET(std::list) expansion = explode(helper1);
+  double res = expansion.minDistSquared(helper2);
   return res;
 }
@@ -137 +204 @@
   res = sqrt(periodicDistanceSquared(point1,point2));
   return res;
+}
+
+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];
 }
 
@@ -145 +248 @@
     M    = new Matrix(*src.M);
     Minv = new Matrix(*src.Minv);
+    conditions = src.conditions;
   }
   return *this;