CWall.h

Go to the documentation of this file.

#ifndef CWALL_H
#define CWALL_H

#include "CParticle.h"
#include <algorithm>
using std::min;



class CWall 
{
 public:
  CWall() {
    indSpecies = 0;
    is_hopper = false;
        isAxisymmetric = false;
        velocity.set_zero();
  }

  void set_is_hopper(){is_hopper=true;}
                
  void set(Vec3D normal_, Mdouble position_) 
  {
    //factor is used to set n to unit length 
    factor = 1. / sqrt(Dot(normal_, normal_));
    normal = normal_ * factor;
    position = position_ * factor;              
    //this makes it a standard wall
    finite_walls.resize(0);
  }

  void is_finite() {finite_walls.resize(0);}
                
  void add_finite_wall(Vec3D normal, Vec3D point) {
          add_finite_wall(normal,Dot(normal,point));
  }

  void add_finite_wall(Vec3D normal_, Mdouble position_) 
  {
    //n is the index of the new wall
    int n = finite_walls.size();
    finite_walls.resize(n+1);
    finite_walls[n].set(normal_,position_);
                
    // AB[n*(n-1)/2+m] is the direction of the intersecting line between walls m and n, m<n
    // A[n*(n-1)/2+m] is a point on the intersecting line between walls m and n, m<n
    // See http://www.netcomuk.co.uk/~jenolive/vect18d.html for finding the line where two planes meet
    AB.resize(n*(n+1)/2);
    A.resize (n*(n+1)/2);
    for(int m=0; m<n; m++) 
      {
        int id = (n-1)*n/2+m;
        //first we cross the wall normals and normalize to obtain AB
        AB[id] = Cross(finite_walls[m].normal, finite_walls[n].normal);
        AB[id] /= sqrt(AB[id].GetLength2());
        //then we find a point A (using AB*x=0 as a third plane)
        Mdouble invdet = 1.0/(+finite_walls[n].normal.X*(finite_walls[m].normal.Y*AB[id].Z-AB[id].Y*finite_walls[m].normal.Z)
                             -finite_walls[n].normal.Y*(finite_walls[m].normal.X*AB[id].Z-finite_walls[m].normal.Z*AB[id].X)
                             +finite_walls[n].normal.Z*(finite_walls[m].normal.X*AB[id].Y-finite_walls[m].normal.Y*AB[id].X));
        A[id] = Vec3D(+(finite_walls[m].normal.Y*AB[id].Z-AB[id].Y*finite_walls[m].normal.Z)*finite_walls[n].position
                      -(finite_walls[n].normal.Y*AB[id].Z-finite_walls[n].normal.Z*AB[id].Y)*finite_walls[m].position
                      +(finite_walls[n].normal.Y*finite_walls[m].normal.Z-finite_walls[n].normal.Z*finite_walls[m].normal.Y)*0.0,
                      -(finite_walls[m].normal.X*AB[id].Z-finite_walls[m].normal.Z*AB[id].X)*finite_walls[n].position
                      +(finite_walls[n].normal.X*AB[id].Z-finite_walls[n].normal.Z*AB[id].X)*finite_walls[m].position
                      -(finite_walls[n].normal.X*finite_walls[m].normal.Z-finite_walls[m].normal.X*finite_walls[n].normal.Z)*0.0,
                      +(finite_walls[m].normal.X*AB[id].Y-AB[id].X*finite_walls[m].normal.Y)*finite_walls[n].position
                      -(finite_walls[n].normal.X*AB[id].Y-AB[id].X*finite_walls[n].normal.Y)*finite_walls[m].position
                      +(finite_walls[n].normal.X*finite_walls[m].normal.Y-finite_walls[m].normal.X*finite_walls[n].normal.Y)*0.0 ) * invdet;
      }
                
    // C[(n-2)*(n-1)*n/6+(m-1)*m/2+l] is a point intersecting walls l, m and n, l<m<n
    C.resize((n-1)*n*(n+1)/6);
    for(int m=0; m<n; m++) 
      {
        for(int l=0; l<m; l++) 
          {
            int id = (n-2)*(n-1)*n/6+(m-1)*m/2+l;
            Mdouble invdet = 1.0/(+finite_walls[n].normal.X*(finite_walls[m].normal.Y*finite_walls[l].normal.Z-finite_walls[l].normal.Y*finite_walls[m].normal.Z)
                                 -finite_walls[n].normal.Y*(finite_walls[m].normal.X*finite_walls[l].normal.Z-finite_walls[m].normal.Z*finite_walls[l].normal.X)
                                 +finite_walls[n].normal.Z*(finite_walls[m].normal.X*finite_walls[l].normal.Y-finite_walls[m].normal.Y*finite_walls[l].normal.X));
            C[id] = Vec3D(+(finite_walls[m].normal.Y*finite_walls[l].normal.Z-finite_walls[l].normal.Y*finite_walls[m].normal.Z)*finite_walls[n].position
                          -(finite_walls[n].normal.Y*finite_walls[l].normal.Z-finite_walls[n].normal.Z*finite_walls[l].normal.Y)*finite_walls[m].position
                          +(finite_walls[n].normal.Y*finite_walls[m].normal.Z-finite_walls[n].normal.Z*finite_walls[m].normal.Y)*finite_walls[l].position,
                          -(finite_walls[m].normal.X*finite_walls[l].normal.Z-finite_walls[m].normal.Z*finite_walls[l].normal.X)*finite_walls[n].position
                          +(finite_walls[n].normal.X*finite_walls[l].normal.Z-finite_walls[n].normal.Z*finite_walls[l].normal.X)*finite_walls[m].position
                          -(finite_walls[n].normal.X*finite_walls[m].normal.Z-finite_walls[m].normal.X*finite_walls[n].normal.Z)*finite_walls[l].position,
                          +(finite_walls[m].normal.X*finite_walls[l].normal.Y-finite_walls[l].normal.X*finite_walls[m].normal.Y)*finite_walls[n].position
                          -(finite_walls[n].normal.X*finite_walls[l].normal.Y-finite_walls[l].normal.X*finite_walls[n].normal.Y)*finite_walls[m].position
                          +(finite_walls[n].normal.X*finite_walls[m].normal.Y-finite_walls[m].normal.X*finite_walls[n].normal.Y)*finite_walls[l].position ) * invdet;
          }
      }
  }
        
  void move(Mdouble position_) {
          position=position_*factor;
          //cout << "pos " << position << endl;
  }

  void move(Vec3D velocity_, Mdouble dt) {
          velocity=velocity_;
          position+=Dot(velocity,normal)*dt;
          //cout << "pos " << position << "vel " << velocity << endl;
  }

  Mdouble get_distance(Vec3D &Position) {return position - Dot(Position, normal);}
                
  bool get_distance_and_normal(Particle &P, Mdouble &distance, Vec3D &normal_return) {
    if (isAxisymmetric) return get_distance_and_normal_for_axissymmetric_wall(P, distance, normal_return);
    static Mdouble distance_new;
    static Mdouble distance2;  
    static Mdouble distance3;
    static int id;  
    static int id2;  
    static int id3;  
                
    if(is_hopper){//this is first becuase i dont know what it does if not
                  // Three posible normals for the hopper, they are averaged
                  // when a particle touch 2 walls
      Vec3D normalOuterCylinder = Vec3D(0,0,0);
      Vec3D normalInnerCylinder = Vec3D(0,0,0);
      Vec3D normalBase = Vec3D(0,0,0);
                  
                  
      //Radial coordinate of the particle
      Mdouble particleRadius = sqrt(sqr(P.get_Position().X) + sqr(P.get_Position().Y));
      Mdouble innerRadius = 20*P.get_Radius();// hardcoded geometry for the moment
      Mdouble outerRadius = 50*P.get_Radius();
      Mdouble distanceOuterCylinder = 0.0;
      Mdouble distanceInnerCylinder = 0.0;
      Mdouble distanceBase = 0.0;
                  
      //cout << "Mass of particles: " << P.get_restitution_coefficient() << endl;
      //set normal outer cilinder
      if(particleRadius + P.get_Radius() > outerRadius){
        normalOuterCylinder = Vec3D(P.get_Position().X, P.get_Position().Y, 0.0)/particleRadius;
        distanceOuterCylinder = outerRadius - particleRadius;
        //cout << "Outer Cylinder "<< endl;
      }

      //\todo set a real hopper with two more walls
     /*  if (P.get_Position().Z - P.Radius < 0.0 && particleRadius +  P.Radius >= innerRadius  */
/*        && P.get_Position().Z > 0){// the particle feels the ground */
/*      normalBase = Vec3D(0.0,0.0,-1.0);//normalized normal outer to the geometry  */
/*      distanceBase = P.get_Position().Z;//the distance to the wall from the center of the particle */
/*      //cout << "Base "<< endl; */
/*       } */
      //\todo set a real hopper with two more walls
      if (abs(P.get_Position().Z) - P.get_Radius() < 0.0 
          && particleRadius > innerRadius){// the particle feels the ground
        
        normalBase = Vec3D(0.0,0.0,-1.0)*SIGN(P.get_Position().Z);
        distanceBase = P.get_Position().Z;
      }
      Vec3D dr = Vec3D(P.get_Position().X, P.get_Position().Y, 0.0)*innerRadius/particleRadius 
        - P.get_Position();
      if(dr.GetLength() < P.get_Radius() ){
        cout << " toy en el borde" << endl;
        distanceBase = dr.GetLength();
        normalBase = dr/distanceBase;
      }
      
      /*   if(particleRadius +  P.Radius > innerRadius  */
      /*         && P.get_Position().Z > -40*P.Radius  */
      /*         && P.get_Position().Z + P.Radius < 0.0 ){//inside wall */
      /*        normalInnerCylinder = Vec3D(-P.get_Position().X, -P.get_Position().Y, 0.0)/particleRadius; */
      /*        distanceInnerCylinder = innerRadius - particleRadius; */
      /*        //cout << "Inner Cylinder "<< endl; */
      /*       } */

      /*  if(P.get_Position().Z > 0.0 && P.get_Position().Z - P.Radius < 0.0 && particleRadius < innerRadius ){// Colliding with line */
      /*        normal_return =  Vec3D(P.get_Position().X,P.get_Position().Y,0)*(innerRadius/particleRadius - 1.0)+ */
      /*          Vec3D(0.0,0.0,-P.get_Position().Z); */
      /*        distance = normal_return.GetLength();  */
      /*        normal_return /= normal_return.GetLength(); */
      /*        if(distance < P.Radius){ */
      /*          cout << "Border normal " << normal_return << " distance " << distance<< endl; */
      /*          return true; */
      /*        } */
      /*       } */

      if(distanceOuterCylinder != 0.0 || distanceBase != 0.0){
        //average normals
                    
        normal_return  = distanceOuterCylinder*normalOuterCylinder 
          + distanceInnerCylinder*normalInnerCylinder 
          + distanceBase*normalBase;
        normal_return /= normal_return.GetLength();
        //  distance = sqrt(SQR(distanceBase) + SQR(distanceOuterCylinder) 
        //          + SQR(distanceInnerCylinder));
        distance = 0.5*(distanceBase + distanceOuterCylinder); 
        //cout << "normal " << normal_return << " distance " << distance<< endl;
        return true;                
      }

      /*     if(distanceBase != 0 || distanceInnerCylinder != 0.0){ //collision with inner walls */
      /*                     //average normals */
                    
      /*                    normal_return  = distanceInnerCylinder*normalInnerCylinder */
      /*                                     + distanceBase*normalBase; */
      /*                    normal_return /= normal_return.GetLength(); */
      /*                    //  distance = sqrt(SQR(distanceBase) + SQR(distanceOuterCylinder) */
      /*                    //      + SQR(distanceInnerCylinder)); */
      /*                    distance = 0.5*(distanceBase + distanceInnerCylinder); */
      /*                    //cout << "normal " << normal_return << " distance " << distance<< endl; */
      /*                    return true; */
      /*                  } */


      //cout << "False " << endl;
      //getchar();
      // no collision, so return false
      return false;               
                  
    }
                
                  
                
    if (!finite_walls.size()) {
      //for standard walls this function sets normal_return and distance and returns true if the particle is in contact
      distance = get_distance(P.get_Position());
      if (distance>=P.get_Radius()) return false;
      normal_return = normal; return true;
    } 

        
                
    //if we are here, this is a finite wall
    distance = -1e20;
    distance2 = -1e20;
    distance3 = -1e20;
    //for all finite walls
    for (unsigned int i=0; i<finite_walls.size(); i++) {
      //calculate the distance to the particle
      distance_new = finite_walls[i].get_distance(P.get_Position());
      //return false if the distance to any one wall is too large (i.e. no contact)
      if (distance_new>=P.get_Radius()) return false;
      //store the minimum distance and the respective wall in "distance" and "id"
      //and store up to two walls (id2, id3) and their distances (distance2, distance3), if the possible contact point is near the intersection between id and id2 (and id3)
      if (distance_new>distance) {
        if (distance>-P.get_Radius()) {
          if (distance2>-P.get_Radius()) {distance3 = distance; id3 = id;}
          else {distance2 = distance; id2 = id;}
        } 
        distance = distance_new; id = i;
      } else if (distance_new>-P.get_Radius()) {
        if (distance2>-P.get_Radius()) {distance3 = distance_new; id3 = i;}
        else {distance2 = distance_new; id2 = i;}
      }
    }
                
    //if we are here, the closest wall is id; 
    //if distance2>-P.Radius (and distance3>-P.Radius), the possible contact point is near the intersection between id and id2 (and id3)
    if (distance2>-P.get_Radius()) {
      //D is the point on wall id closest to P
      Vec3D D = P.get_Position() + finite_walls[id].normal * distance;
      //If the distance of D to id2 is positive, the contact is with the intersection
      bool intersection_with_id2 = (finite_walls[id2].get_distance(D)>0.0);

      if (distance3>-P.get_Radius()&&(finite_walls[id3].get_distance(D)>0.0)) {
        if (intersection_with_id2) { //possible contact is with intersection of id,id2,id3
          //we know id2<id3
          int index = 
            (id<id2)?( (id3-2)*(id3-1)*id3/6+(id2-1)*id2/2+id  ):
            (id<id3)?( (id3-2)*(id3-1)*id3/6+(id -1)*id /2+id2 ):
            ( (id -2)*(id -1)*id /6+(id3-1)*id3/2+id2 );
          normal_return = P.get_Position() - C[index];
          distance = sqrt(normal_return.GetLength2());
          if (distance>=P.get_Radius()) return false; //no contact
          normal_return /= -distance; 
          return true; //contact with id,id2,id3
        } else { intersection_with_id2 = true; distance2 = distance3; id2 = id3; }
      }
                        
      if (intersection_with_id2) { //possible contact is with intersection of id,id2
        int index = (id>id2)?((id-1)*id/2+id2):((id2-1)*id2/2+id);
        Vec3D AC = P.get_Position() - A[index];
        normal_return = AC - AB[index] * Dot(AC,AB[index]);
        distance = sqrt(normal_return.GetLength2());
        if (distance>=P.get_Radius()) return false; //no contact
        normal_return /= -distance; 
        return true; //contact with two walls
      }
    }
    //contact is with id
    normal_return = finite_walls[id].normal; 
    return true;
                
  }

  friend inline std::ostream& operator<<(std::ostream& os, CWall &w) 
  {
    os << "numFiniteWalls " << w.finite_walls.size(); 
    if (w.finite_walls.size()) {
        for (vector<CWall>::iterator it = w.finite_walls.begin(); it!=w.finite_walls.end(); ++it) {
          os << " normal " << it->normal << " position " << it->position;
        }
    } else {
        os << " normal " << w.normal << " position " << w.position;
        //optional output
        if (w.velocity.GetLength2()) os << " velocity " << w.velocity;
    }
    return os;
  }
        
  friend inline std::istream& operator>>(std::istream& is, CWall &w) 
  {
    int n;
    is >> n;
    if (n!=0) {
        Vec3D normal;
        Mdouble position;
        for (int i=0; i<n; i++) {
          is >> normal >> position;
          w.add_finite_wall(normal,position);
        }               
    } else {
        is >> w.normal >> w.position;
    }
    return is;
  }

    void read(std::istream& is)  { 
        string dummy;
        int n;
        is >> dummy >> n;
        if (n!=0) {
            Vec3D normal;
            Mdouble position;
            for (int i=0; i<n; i++) {
                is >> dummy >> normal >> dummy >> position;
                add_finite_wall(normal,position);
            }           
        } else {
            //use stringstream to ensure that old restart files w/o velocity don't break the code
            stringstream ss("");
            string s;
            getline(is,s);
            ss << s;
            ss >> dummy >> normal >> dummy >> position >> dummy >> velocity;

        }
    }
        
    void print(std::ostream& os) { 
        if (finite_walls.size()) {
        os << "finite_wall( "; 
        for (vector<CWall>::iterator it = finite_walls.begin(); it!=finite_walls.end(); ++it)
            os << "normal:" << it->get_normal() << ", position:" << it->get_position() << " "; 
            os << " )"; 
        } else os << "wall( normal:" << normal << ", position:" << position << ")"; 
    }
        
    Vec3D get_normal() {return normal;}
    void set_normal(Vec3D new_) {normal=new_;}

    Mdouble get_position() {return position;}


    Vec3D get_velocity() {return velocity;}
    void set_velocity(Vec3D new_) {velocity = new_;}

    bool get_distance_and_normal_for_axissymmetric_wall(Particle &P, Mdouble &distance, Vec3D &normal_return) {
            //transform to axisymmetric coordinates
            Vec3D PO = P.get_Position() - axisOrigin;
            P.get_Position().Z = Dot(axisOrientation,PO);
            P.get_Position().Y = 0.0;
            Vec3D tangential = PO-P.get_Position().Z*axisOrientation;
            P.get_Position().X = tangential.GetLength();
            tangential /= P.get_Position().X;
            Vec3D normal_axisymmetric_coordinates;
            //determine wall distance, normal and contact in axosymmetric coordinates
            isAxisymmetric = false;
            bool retval = get_distance_and_normal(P, distance, normal_axisymmetric_coordinates);
            isAxisymmetric = true;
            //transform from axisymmetric coordinates
            normal_return = normal_axisymmetric_coordinates.Z * axisOrientation 
                    + tangential*normal_axisymmetric_coordinates.X;             
            P.set_Position(PO + axisOrigin);
            return retval;
    }

    void setSymmetryAxis(Vec3D new_axisOrigin, Vec3D new_axisOrientation) {
            axisOrigin = new_axisOrigin;
            axisOrientation = new_axisOrientation;
            axisOrientation.normalize();
            isAxisymmetric = true;
            cout 
            << "axisOrigin" << axisOrigin
            << "axisOrientation" << axisOrientation
            << "isAxisymmetric" << isAxisymmetric
            << endl;
    } 
    
    void removeSymmetryAxis() {
            isAxisymmetric = false;
    }           
                
private:
  vector<CWall> finite_walls;
  vector<Vec3D> A;  //<A[n*(n-1)/2+m] is a point on the intersecting line between walls m and n, m<n
  vector<Vec3D> AB; //<AB[n*(n-1)/2+m] is the direction of the intersecting line between walls m and n, m<n
  vector<Vec3D> C;  //<C[(n-2)*(n-1)*n/6+(m-1)*m/2+l] is a point intersecting walls l, m and n, l<m<n
        
  Vec3D normal;     //<outward unit normal vector
  Mdouble position;  //<position n*x=p
  Mdouble factor;    //<This is the normal to rescale to unit vectoers.


 public:
  Vec3D velocity;        
  int indSpecies;
  bool is_hopper;
 
 
  bool isAxisymmetric;
  Vec3D axisOrigin;
  Vec3D axisOrientation;

};


class CWallPeriodic 
{
public:

        CWallPeriodic() {}

        CWallPeriodic(const CWallPeriodic &p) {
                left_wall = p.left_wall;
                normal = p.normal;
                position_left = p.position_left;
                position_right = p.position_right;
                factor = p.factor;
                shift = p.shift;
        }

  void set(Vec3D normal_, Mdouble position_left_, Mdouble position_right_) {
    // factor is used to set normal to unit length
    factor = 1. / sqrt(Dot(normal_, normal_));
    normal = normal_ * factor;
    position_left = position_left_ * factor;
    position_right = position_right_ * factor;
    shift = normal * (position_right - position_left);
  }

  void set(Vec3D normal_, Mdouble position_left_, Mdouble position_right_, Vec3D shift_direction) {
    // factor is used to set normal to unit length
    factor = 1. / sqrt(Dot(normal_, normal_));
    normal = normal_ * factor;
    position_left = position_left_ * factor;
    position_right = position_right_ * factor;
    // factor is used to set shift vector to correct length
    factor = (position_right - position_left) * Dot(shift_direction, normal);
    shift = shift_direction * factor;
  }
  
  void move_left(Mdouble position_)
  {
          position_left=position_*factor;
          shift = normal * (position_right - position_left);
  }
  
  void move_right(Mdouble position_)
  {
          position_right=position_*factor;
          shift = normal * (position_right - position_left);
  }


  Mdouble distance(Particle &P) {
    return distance(P.get_Position());
  }

  Mdouble distance(Vec3D &P) {
    Mdouble position = Dot(P, normal);

    if (position - position_left < position_right - position) {
      left_wall = true;
      return position - position_left;
    } else {
      left_wall = false;
      return position_right - position;
    }
  }

        void shift_position(Vec3D &Position) {
                if (left_wall) {
                        Position += shift;
                        left_wall = false;
                }
                else {
                        Position -= shift;
                        left_wall = true;
                }
        }

        Vec3D get_shifted_position(Vec3D &Position) {
                if (left_wall) {
                        return Position + shift;
                }
                else {
                        return Position - shift;
                }
        }
        
        void shift_positions(Vec3D &PI, Vec3D &PJ) {
                if (left_wall) {
                        PI += shift;
                        PJ += shift;
                        left_wall = false;
                }
                else {
                        PI -= shift;
                        PJ -= shift;
                        left_wall = true;
                }
        }
                
        void get_close_together(Vec3D &P,Vec3D &Q) {
                Mdouble PQdotn = Dot(P-Q, normal);
                Mdouble shift_norm2 = shift.GetLength2();
                //Check if P is so far from Q that a shift would move it closer
                if (sqr(PQdotn) > .25 * shift_norm2) {
                        //Now determine the direction of the shift
                        if (PQdotn>0.0) P -= shift;
                        else P += shift;
                }
        }
                
        friend inline std::ostream& operator<<(std::ostream& os, const CWallPeriodic &w) 
        {
                os<< "normal " << w.normal 
      << " pos_left " << w.position_left 
      << " pos_right " << w.position_right 
      << " shift " << w.shift;
                return os;
        }
        
        friend inline std::istream& operator>>(std::istream& is, CWallPeriodic &w) 
        {
                is >> w.normal >> w.position_left >> w.position_right >> w.shift;
                return is;
        }

        void read(std::istream& is)  { 
    string dummy;
                is >> dummy >> normal >> dummy >> position_left >> dummy >> position_right >> dummy >> shift;
        }
        
        void print(std::ostream& os, bool print_all=false) {
                os << "periodic_wall( normal:" << normal 
                         << ", position_left:" << position_left 
                         << ", position_right:" << position_right;
                if (print_all) os       << ", shift:" << shift << ", left_wall:" << left_wall;
                os << ")";
        }
                
        Vec3D& get_normal() {return normal;}
        void set_normal(Vec3D new_) {normal=new_;}
        //~ Vec3D get_normal() const {return normal;}
        bool get_left_wall() const {return left_wall;}

 private:
  bool left_wall;        
  Vec3D normal;          
  Mdouble position_left;  
  Mdouble position_right; 
  Mdouble factor;        
  Vec3D shift;           
};

#endif