ESyS-Particle
4.0.1
|
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