/Users/jack/Code/basso_dev/inc/ElasticTruss3D.h

Go to the documentation of this file.

#ifndef _ELASTIC_TRUSS_3D_BASSO_H_
#define _ELASTIC_TRUSS_3D_BASSO_H_

// std includes

// Basso includes
#include "basso.h"
#include "StructuralFormulation.h"
#include "MatElastic.h"
#include "basic_fem_operations.h"

namespace basso {



class ElasticTruss3D : public StructuralFormulation {

public:
        
        // PUBLIC DEFS
        
        // CONSTRUCTORS 
        ElasticTruss3D( list< Element > &elem, const Array< Node > &nodes, DofMap &dofmap, 
                const Material *material, const Numeric a ) 
                : StructuralFormulation( elem, nodes, dofmap, kSTRESS1D )
        { 
                this->mSdim = 3;
                this->fMat = material; 
                Initialize(a);
                CheckElements();
        }
        
        ~ElasticTruss3D() { }
        
        // ACCESSORS
        
        // MEMBER FUNCTIONS

        
protected:
        virtual void GetElementStiffnessMatrix( nMatrix &ke, list<Element>::const_iterator eItr  );
        
        virtual void GetElementMassMatrix( nMatrix &me, list<Element>::const_iterator eItr  );
        
        void CheckElements() const;
        
        virtual void Initialize( const Numeric a );
        
private:
        Numeric mArea;

};

void ElasticTruss3D::GetElementStiffnessMatrix( nMatrix &ke, list<Element>::const_iterator eItr  )  
{
        Numeric x1=this->fNode[ eItr->Node(0) ].x(), x2=this->fNode[ eItr->Node(1) ].x(),
                        y1=this->fNode[ eItr->Node(0) ].y(), y2=this->fNode[ eItr->Node(1) ].y(), 
                        z1=this->fNode[ eItr->Node(0) ].z(), z2=this->fNode[ eItr->Node(1) ].z();
        Numeric L=sqrt( pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2) );
        Numeric k=mArea*fMat->YoungsModulus()/L;

        Numeric ls=(x2-x1)/L, ms=(y2-y1)/L, ns=(z2-z1)/L;
        ke(0,0)=k*ls*ls; ke(0,1)=k*ls*ms; ke(0,2)=k*ls*ns; ke(0,3)=-k(0,0); ke(0,4)=-ke(0,1); ke(0,5)=-ke(0,2);
        ke(1,0)=ke(0,1); ke(1,1)=k*ms*ms; ke(1,2)=k*ms*ns; ke(1,3)=-k(1,0); ke(1,4)=-ke(1,1); ke(1,5)=-ke(1,2);
        ke(2,0)=ke(0,2); ke(2,1)=ke(1,2); ke(2,2)=k*ns*ns; ke(2,3)=-k(2,0); ke(2,4)=-ke(2,1); ke(2,5)=-ke(2,2);
        ke(3,0)=ke(0,3); ke(3,1)=ke(1,3); ke(3,2)=ke(2,3); ke(3,3)= k(0,0); ke(3,4)= ke(0,1); ke(3,5)= ke(0,2);
        ke(4,0)=ke(0,4); ke(4,1)=ke(1,4); ke(4,2)=ke(2,4); ke(4,3)= k(3,4); ke(4,4)= ke(1,1); ke(4,5)= ke(1,2);
        ke(5,0)=ke(0,5); ke(5,1)=ke(1,5); ke(5,2)=ke(2,5); ke(5,3)= k(3,5); ke(5,4)= ke(4,5); ke(5,5)= ke(2,2);
        
}

void ElasticTruss3D::GetElementMassMatrix( nMatrix &me, list<Element>::const_iterator eItr  )  
{
        Numeric L=sqrt( pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2) );
        Numeric mn=0.5*mArea*L*fMat->Density();

        me(0,0)=mn/3;    me(0,1)=0.0;     me(0,2)=0.0;     me(0,3)=mn/6;    me(0,4)=0.0;     me(0,5)=0.0; 
        me(1,0)=0.0;     me(1,1)=me(0,0); me(1,2)=0.0;     me(1,3)=0.0;     me(1,4)=me(0,2); me(1,5)=0.0;
        me(2,0)=0.0;     me(2,1)=0.0;     me(2,2)=me(0,0); me(2,3)=0.0;     me(2,4)=0.0;     me(2,5)=me(0,2); 
        me(3,0)=me(0,2); me(3,1)=0.0;     me(3,2)=0.0;     me(3,3)=me(0,0); me(3,4)=0.0;     me(3,5)=0.0; 
        me(4,0)=0.0;     me(4,1)=me(0,2); me(4,2)=0.0;     me(4,3)=0.0;     me(4,4)=me(0,0); me(4,5)=0.0;
        me(5,0)=0.0;     me(5,1)=0.0;     me(5,2)=me(0,2); me(5,3)=0.0;     me(5,4)=0.0;     me(5,5)=me(0,0);
        
}

void ElasticTruss3D::Initialize( const Numeric a ) 
{  
        this->mArea = a; 
        resize( this->activeLocalDofs, this->mSdim );
        this->activeLocalDofs[0]=kDISPx; 
        this->activeLocalDofs[1]=kDISPy;
        this->activeLocalDofs[2]=kDISPz;
}

void ElasticTruss3D::CheckElements() const
{
        list< Element >::const_iterator eItr;
        for ( eItr=this->fElement->begin(); eItr!=this->fElement->end(); ++eItr ) 
                if ( eItr->Type()!=kLINE2 )
                        WarningMessage("ElasticTruss3D","Elements must be kLINE2; other elements will be ignored");
}



} // end namespace

#endif

---