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

Go to the documentation of this file.

#ifndef _STRUCTURAL_FORUMULATION_BASSO_H_
#define _STRUCTURAL_FORUMULATION_BASSO_H_

// std includes
#include <list>

// Basso includes
#include "basso.h"
#include "Formulation.h"
#include "StressState.h"

namespace basso {

        enum StressStateType { kSTRESS_1D, kPLANE_STRAIN, kPLANE_STRESS, kAXISYMMETRIC, kSTRESS_3D };

class StructuralFormulation : public Formulation {

public:
        
        // PUBLIC DEFS
        
        // CONSTRUCTORS 
        StructuralFormulation( list< Element > &elem, const Array< Node > &nodes, StressStateType ss=kSTRESS_3D ) 
                : Formulation( elem, nodes ) 
                { 
                        switch (ss)
                        {
                                case kSTRESS_1D:
                                fStressState = new StressState1D;
                                break;
                                
                                case kPLANE_STRAIN:
                                fStressState = new PlaneStrainState;
                                break;
                                
                                case kPLANE_STRESS:
                                fStressState = new PlaneStressState;
                                break;
                                
                                case kAXISYMMETRIC:
                                fStressState = new AxisymmetricState;
                                break;
                                
                                default:
                                fStressState = new StressState3D;
                                break;
                                
                        }
                        this->mVdim = this->fStressState->VoigtDimension();
                        resize( this->mCmat, this->mVdim, this->mVdim ); 
                }
        
        virtual ~StructuralFormulation() { }
        
        // ACCESSORS
        
        // MEMBER FUNCTIONS
                
        virtual int NumLocalDofs() const { return this->mSdim; }

        virtual void NodalAvgStress( const nArray &d, Array< nArray > &stress, Array< nArray > &strain );

protected:
        virtual void SetCMatrix() = 0;
                
        virtual void GetElementStiffnessMatrix( nMatrix &ke, list< Element >::const_iterator eItr  );

        virtual Numeric SetBMatrix( list< Element >::const_iterator eItr, const Point &pt )
        {
                Numeric jac = this->SetDNx( eItr, pt );
                form_bmatrix( this->mDNx, this->mBmat );
                return jac;             
        }

protected:
        
        virtual bool SetElementParameters( list<Element>::const_iterator eItr, int qorder=2 );
        
protected:
        
        const Material *fMat;  
        StressState *fStressState;
        
        nMatrix mBmat, mCmat;
        int mVdim, mEdim;

};

void StructuralFormulation::NodalAvgStress( 
        const nArray &d, Array< nArray > &stress, Array< nArray > &strain )  
{
        if  ( fDofMap==0 )
        {
                WarningMessage("StructuralFormulation::NodalAvgStress","DofMap not set");
                return;
        }
        
        if ( stress.size()!=d.size() )
                resize( stress, d.size() );
        if ( strain.size()!=d.size() )
                resize( strain, d.size() );
        
        iArray navg( d.size() ); clear( navg );
        nArray de, estrain(mVdim), estress(mVdim);
        Array<Point> pts;
        list< Element >::const_iterator eItr;
        
        for ( eItr=this->fElement->begin(); eItr!=this->fElement->end(); ++eItr )  // element loop
        {
        
                if ( SetElementParameters( eItr ) ) {
                        resize( de, mNne*NumLocalDofs() );
                        eItr->ParentCoord( pts );
                }
                
                fDofMap->Gather( eItr->Connectivity(), ActiveDofs(), d, de );
                for ( int n=0; n<this->mNne; ++n )
                {               
                        SetBMatrix( eItr, pts[n] );
                        mult( mBmat, de, estrain );
                        mult( mCmat, estrain, estress );
                        int I=eItr->NodeId(n);
                        if ( navg[ I ]==0 ) // put stress and strain in 
                        {
                                if ( stress[I].size()!=estress.size() )
                                        resize( stress[I], estress.size() );
                                if ( strain[I].size()!=estrain.size() )
                                        resize( strain[I], estrain.size() );
                                stress[I]=estress;
                                strain[I]=estrain;
                        }
                        else
                        {
                                stress[I]=estress;
                                strain[I]=estrain;
                        }
                        ++navg[I];
                        
                }
        }
        
        for ( int I=0; I<d.size(); ++I )
                for ( int i=0; i<stress[I].size(); ++i )
                {
                        stress[I][i] = stress[I][i]/navg[I];
                        strain[I][i] = strain[I][i]/navg[I];
                }
                        
}

bool StructuralFormulation::SetElementParameters( list<Element>::const_iterator eItr, int qorder )
{
        if ( this->mLastElemType == eItr->Type() )  // do nothin
                return false;
        this->mNne = eItr->NumNodes();
        this->mEdim = eItr->Dimension();
        resize( this->mDNxi, this->mNne, this->mEdim );
        resize( this->mDNx, this->mNne, this->mSdim ); 
        resize( this->mNa, this->mNne );
        resize( this->mCoord, this->mNne, this->mSdim );  
        resize( this->mBmat, this->mVdim, this->mSdim*this->mNne );
        eItr->Quadrature( qorder, this->mQuadrature );
        this->mLastElemType = eItr->Type();
        return true;
}

void StructuralFormulation::GetElementStiffnessMatrix( nMatrix &ke, list<Element>::const_iterator eItr  )  
{
        SetElementParameters( eItr, 2*(eItr->Order()-1) );
        int kdim = eItr->NumNodes()*NumLocalDofs();
        resize( ke, kdim, kdim );
        clear(ke);
        Numeric jac;
        QuadratureRule::const_iterator qItr;
        SetCMatrix();
        for ( qItr=QuadratureBegin(); qItr!=QuadratureEnd(); ++qItr )   // quadrature loop
        {
                jac = SetBMatrix( eItr, qItr->Pt() );
                femmult_add( transposed(this->mBmat), mCmat, this->mBmat, jac*qItr->Wt(), ke );
        }               
        
}

} // end namespace

#endif

---