ESyS-Particle  4.0.1
SoftBWallInteraction.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 
00014 #ifndef MODEL_SOFTBWALLINTERACTION_HPP
00015 #define MODEL_SOFTBWALLINTERACTION_HPP
00016 
00017 #include "SoftBWallInteraction.h"
00018 
00019 template <class T>
00020 CSoftBondedWallInteraction<T>::CSoftBondedWallInteraction(T* p,CWall* w,double normalK,double shearK,bool scaling,bool iflag):
00021   AWallInteraction<T>(p,w,iflag)
00022 {
00023   double scale;
00024 
00025   if (scaling) {
00026     // scale stiffness to particle cross section
00027     if(CParticle::getDo2dCalculations()){ // 2D
00028 //      scale=2.0*this->m_p->getRad();
00029       scale=1.0;
00030     } else { // 3D
00031 //      scale=3.1415926536*this->m_p->getRad()*this->m_p->getRad();
00032       scale=3.1415926536*this->m_p->getRad();
00033     }
00034   }
00035   else {
00036      scale = 1.0;
00037   }
00038 
00039   m_normalK=normalK*scale;
00040   m_shearK=shearK*scale;
00041 }
00042 
00046 template <class T>
00047 void CSoftBondedWallInteraction<T>::calcForces()
00048 {
00049   const Vec3 &n  = this->m_wall->getNormal();
00050   const Vec3 relDisplacement =
00051     (
00052       this->m_p->getTotalDisplacement()
00053       -
00054       this->m_wall->getTotalDisplacement()
00055     );
00056   
00057   const double normalDist = dot(relDisplacement, n)/(n.norm());
00058   const Vec3 normalForce = ((m_normalK*normalDist)/n.norm())*n;
00059 
00060   /*
00061    * Tangent displacement vector is relDisplacement point
00062    * projected onto plane.
00063    */
00064   const Vec3 tangentDisplacement =
00065     (relDisplacement - ((normalDist/(n.norm()))*n));
00066   const Vec3 tangentForce = m_shearK*tangentDisplacement;
00067 
00068   const Vec3 totalForce = normalForce + tangentForce;
00069   this->m_p->applyForce(-1.0*totalForce,this->m_p->getPos());
00070   if(this->m_inner_flag) this->m_wall->addForce(totalForce);
00071 }
00072 
00073 #endif