AngledPeriodicBoundary.h

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef AngledPeriodicBoundary_H
#define AngledPeriodicBoundary_H

#include "../BaseBoundary.h"
#include "../Particles/TangentialSpringParticle.h"
#include "../Particles/BaseParticle.h"

class AngledPeriodicBoundary : public BaseBoundary {
public:

    virtual AngledPeriodicBoundary* copy() const
    {
        #ifdef CONSTUCTOR_OUTPUT
            std::cerr << "virtual AngledPeriodicBoundary* copy() const finished" << std::endl;
        #endif                
        return new AngledPeriodicBoundary(*this);
    }
    
    
    //todo constructors instead of set functions?
    //AngledPeriodicBoundary (Vec3D normal_left_, Vec3D normal_right_, Vec3D origin_) 

    void set (Vec3D normal_left_, Vec3D normal_right_, Vec3D origin_) {
        origin=origin_;
        normal_left = normal_left_ / GetLength(normal_left_);
        normal_right = normal_right_ / GetLength(normal_right_);
        common_axis = Cross(normal_left, normal_right);
        common_axis /= GetLength(common_axis);
        radialAxis_left = Cross(normal_left, common_axis);
        radialAxis_right = Cross(normal_right, common_axis);
        diff_normal = normal_right-normal_left;
        diff_radial = radialAxis_right-radialAxis_left;
        //angularShift = 0;
        // std::cout << "common_axis " << common_axis
        // << ", radialAxis_left " << radialAxis_left
        // << ", radialAxis_right " << radialAxis_right
        // << ", angularShift " << angularShift
        // << std::endl; 
    }
    Mdouble distance(BaseParticle &P) {
        return distance(P.get_Position());
    }

    //this function should be cheap, as it has to be computed for all particles
    Mdouble distance(const Vec3D &P) {
        Vec3D position = P-origin;
        Mdouble distance_left = Dot(position,normal_left);
        Mdouble distance_right = -Dot(position,normal_right);

        if (distance_left<distance_right) {
            left_wall = true;
            //std::cout << "left wall, " << position << ", distance " << distance_left << "<" << distance_right << std::endl;
            return distance_left;
        } else {
            left_wall = false;
            //std::cout << "right wall, " << position << ", distance " << distance_right << "<" << distance_left << std::endl;
            return distance_right;
        }
    }

    void shift_position(BaseParticle* P) {
        Vec3D position = P->get_Position()-origin;
        if (left_wall) {
            Mdouble normalDistance = Dot(position,normal_left);
            Mdouble radialDistance = Dot(position,radialAxis_left);
            P->move(normalDistance*diff_normal+radialDistance*diff_radial);
            Mdouble normalVelocity = Dot(P->get_Velocity(),normal_left);
            Mdouble radialVelocity = Dot(P->get_Velocity(),radialAxis_left);
            P->accelerate(normalVelocity*diff_normal+radialVelocity*diff_radial);
            Mdouble normalAngularDistance = Dot(P->get_Angle(),normal_left);
            Mdouble radialAngularDistance = Dot(P->get_Angle(),radialAxis_left);
            P->rotate(normalAngularDistance*diff_normal+radialAngularDistance*diff_radial);
            left_wall = false;
            //std::cout << "shift to right wall, " << P->get_Position() << std::endl;
        }
        else {
            Mdouble normalDistance = Dot(position,normal_right);
            Mdouble radialDistance = Dot(position,radialAxis_right);
            P->move(-normalDistance*diff_normal-radialDistance*diff_radial);
            Mdouble normalVelocity = Dot(P->get_Velocity(),normal_right);
            Mdouble radialVelocity = Dot(P->get_Velocity(),radialAxis_right);
            P->accelerate(-normalVelocity*diff_normal-radialVelocity*diff_radial);
            Mdouble normalAngularDistance = Dot(P->get_Angle(),normal_right);
            Mdouble radialAngularDistance = Dot(P->get_Angle(),radialAxis_right);
            P->rotate(-normalAngularDistance*diff_normal-radialAngularDistance*diff_radial);
            left_wall = true;
            //std::cout << "shift to left wall, " << P->get_Position() << std::endl;
        }
    }

    // ///returns the shifted particle w/o actually shifting
    // Vec3D get_shifted_position(Vec3D &Position) {
    //     if (left_wall) {
    //         return Position + shift;
    //     }
    //     else {
    //         return Position - shift;
    //     }
    // }
    
    // ///shifts two particles
    // 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;
    //     }
    // }
        
    // ///shift P such that it is closest to Q
    // 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;
    //     }
    // }
        

    void read(std::istream& is) { 
        std::string dummy;
        is >> dummy >> normal_left >> dummy >> normal_right >> dummy >> origin;
    }
    
    void print(std::ostream& os) const {
        os << "AngledPeriodicBoundary normal_left " << normal_left 
            << " normal_right " << normal_right 
            << " origin " << origin;
    }
        
    Vec3D& get_normal() {
        if (left_wall) return normal_left;
        else return normal_right;
    }
    
    int createPeriodicParticles(BaseParticle *P, ParticleHandler &pH)
    {
        //std::cout << "createPeriodicParticles" << std::endl;
        if (distance(*P)<P->get_InteractionRadius()+pH.getLargestParticle()->get_InteractionRadius())
        {
            BaseParticle* F0=P->copy();
            shift_position(F0);
                        
            //If the Particle includes TangentalSprings reverse them
            TangentialSpringParticle* TSParticle=dynamic_cast<TangentialSpringParticle*>(F0);
            if(TSParticle) {
                TSParticle->reverseTangentialSprings();
                for(std::vector<CTangentialSpring>::iterator it = TSParticle->get_TangentialSprings().begin(); it!=TSParticle->get_TangentialSprings().end();it++)
                {
                    //std::cout << it->delta << std::endl;
                    if (!left_wall) {
                        Mdouble normalDistance = Dot(it->delta,normal_left);
                        Mdouble radialDistance = Dot(it->delta,radialAxis_left);
                        it->delta += normalDistance*diff_normal+radialDistance*diff_radial;
                        normalDistance = Dot(it->RollingSpring,normal_left);
                        radialDistance = Dot(it->RollingSpring,radialAxis_left);
                        it->RollingSpring += normalDistance*diff_normal+radialDistance*diff_radial;
                        normalDistance = Dot(it->TorsionSpring,normal_left);
                        radialDistance = Dot(it->TorsionSpring,radialAxis_left);
                        it->TorsionSpring += normalDistance*diff_normal+radialDistance*diff_radial;
                    } else {
                        Mdouble normalDistance = Dot(it->delta,normal_right);
                        Mdouble radialDistance = Dot(it->delta,radialAxis_right);
                        it->delta -= normalDistance*diff_normal+radialDistance*diff_radial;
                        normalDistance = Dot(it->RollingSpring,normal_right);
                        radialDistance = Dot(it->RollingSpring,radialAxis_right);
                        it->RollingSpring -= normalDistance*diff_normal+radialDistance*diff_radial;
                        normalDistance = Dot(it->TorsionSpring,normal_right);
                        radialDistance = Dot(it->TorsionSpring,radialAxis_right);
                        it->TorsionSpring -= normalDistance*diff_normal+radialDistance*diff_radial;
                    }
                    //std::cout << it->delta << std::endl;
                }        
            }
    
            //If Particle is Mdouble shifted, get correct original particle
            BaseParticle* From=P;
            while(From->get_PeriodicFromParticle()!=NULL)
                From=From->get_PeriodicFromParticle();        
            F0->set_periodicFromParticle(From);

            pH.addObject(F0);
            return 1;
        }
        return 0;
    }        

    bool checkBoundaryAfterParticleMoved(BaseParticle *P, ParticleHandler &pH UNUSED)
    {
        //std::cout << "checkBoundaryAfterParticleMoved" << std::endl;
        if (distance(*P)<0)
            shift_position(P);
        return false;
    }
    
 private:
    //values set by the user
    Vec3D normal_left;         
    Vec3D normal_right;        
    Vec3D origin;                    
    //values set by the code
    bool left_wall;                
    Vec3D radialAxis_left;         
    Vec3D radialAxis_right;        
    Vec3D diff_radial;
    Vec3D diff_normal;   
    Vec3D common_axis;
    Vec3D angularShift;
};
#endif