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

Go to the documentation of this file.

#ifndef _FORUMULATION_BASSO_H_
#define _FORUMULATION_BASSO_H_

// std includes
#include <list>

// Basso includes
#include "basso.h"
#include "Element.h"
#include "Node.h"
#include "Material.h"
#include "DofMap.h"
#include "basic_fem_operations.h"

namespace basso {


class Formulation {

public:
        
        // PUBLIC DEFS
        
        // CONSTRUCTORS 
        Formulation( list< Element > &elem, const Array< Node > &nodes ) 
        { 
                fDofMap = 0;
                fElement = &elem; fNode = &nodes;
                mLastElemType = kNONE; 
                Initialize();
        }
        
        virtual ~Formulation() {  }
        
        // ACCESSORS
        int NumNodes() const { return vect_size(*fNode); }
        
        int NumElements() const { return fElement->size(); }
        
        virtual int NumLocalDofs() const = 0;
        
        virtual void ActivateDofs( DofMap &dofmap );
        
        virtual const Array<Dof> &ActiveDofs() const { return activeLocalDofs; }
                
        virtual void ContributeStiffnessMatrix( DynamicCSRMatrix &K );
        
        virtual void ContributeMassMatrix( DynamicCSRMatrix &M );
        
protected:
        virtual void GetElementStiffnessMatrix( nMatrix &ke, list<Element>::const_iterator eItr  );
        
        virtual void GetElementMassMatrix( nMatrix &me, list<Element>::const_iterator eItr  );

        QuadratureRule::const_iterator QuadratureBegin() const { return mQuadrature.begin(); }
        QuadratureRule::const_iterator QuadratureEnd() const { return mQuadrature.end(); }

        list<Element>::const_iterator ElementBegin() const { return fElement->begin(); }
        list<Element>::const_iterator ElementEnd() const { return fElement->end(); }

        virtual Numeric SetDNx( list<Element>::const_iterator eItr, const Point &pt );
        
        void Initialize(  ) { }
        
        virtual bool SetElementParameters( list<Element>::const_iterator eItr, int qorder=2 );
        
        virtual Numeric MassFactor() const { return 1.0; }
        
        virtual Numeric StiffnessFactor() const { return 1.0; }

protected:
        list< Element > *fElement;
        const Array< Node > *fNode;
        Array<Dof> activeLocalDofs;
        
        nMatrix mDNxi, mDNx, mCoord;
        nArray mNa;
        int mNne, mSdim;
        QuadratureRule mQuadrature;
        
        DofMap *fDofMap;
        
        BasisType mLastElemType;

};

void Formulation::ActivateDofs( DofMap &dofmap )  
{
        fDofMap=&dofmap;
        list< Element >::const_iterator eItr;
        Element::const_iterator nItr;
        for ( eItr=this->fElement->begin(); eItr!=this->fElement->end(); ++eItr ) 
                for ( nItr=eItr->begin(); nItr!=eItr->end(); ++nItr ) 
                        for ( int s=0; s<NumLocalDofs(); ++s )
                                fDofMap->AddNodalDof( *nItr, this->activeLocalDofs[s] );
        fDofMap->RenumberGlobalDofs();          
}

void Formulation::ContributeMassMatrix( DynamicCSRMatrix &M )  
{
        if  ( fDofMap==0 )
        {
                WarningMessage("Formulation::ContributeMassMatrix","DofMap not set");
                return;
        }
        nMatrix me;

        list< Element >::const_iterator eItr;
        QuadratureRule::const_iterator qItr;
        Element::const_iterator nItr;
        for ( eItr=this->fElement->begin(); eItr!=this->fElement->end(); ++eItr )  // element loop
        {
                GetElementMassMatrix( me, eItr );
                this->fDofMap->Scatter( eItr->Connectivity(), activeLocalDofs, me, M );
        }
}

void Formulation::ContributeStiffnessMatrix( DynamicCSRMatrix &K )  
{
        if  ( fDofMap==0 )
        {
                WarningMessage("Formulation::ContributeStiffnessMatrix","DofMap not set");
                return;
        }
        nMatrix ke;

        list< Element >::const_iterator eItr;
        QuadratureRule::const_iterator qItr;
        Element::const_iterator nItr;
        for ( eItr=this->fElement->begin(); eItr!=this->fElement->end(); ++eItr )  // element loop
        {
                GetElementStiffnessMatrix( ke, eItr );
                this->fDofMap->Scatter( eItr->Connectivity(), activeLocalDofs, ke, K );
        }               
}

void Formulation::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;
        for ( qItr=QuadratureBegin(); qItr!=QuadratureEnd(); ++qItr )   // quadrature loop
        {
                jac = SetDNx( eItr, qItr->Pt() );
                femmult_add( this->mDNx, transposed( this->mDNx ), StiffnessFactor()*jac*qItr->Wt(), ke );
        }               
        
}

void Formulation::GetElementMassMatrix( nMatrix &me, list<Element>::const_iterator eItr  )  
{
        SetElementParameters( eItr, 2*(eItr->Order()) );
        int mdim = eItr->NumNodes()*NumLocalDofs();
        resize( me, mdim, mdim );
        clear(me);
        Numeric jac;
        QuadratureRule::const_iterator qItr;
        for ( qItr=QuadratureBegin(); qItr!=QuadratureEnd(); ++qItr )   // quadrature loop
        {       
                eItr->Na( qItr->Pt(), this->mNa );
                eItr->DNa( qItr->Pt(), this->mDNxi );
                jac = jacobian( this->mCoord, this->mDNxi );
                for ( int I=0; I<mNne; ++I )
                        for ( int J=0; J<mNne; ++J )
                                for ( int s=0; s<NumLocalDofs(); ++s )
                                        me( NumLocalDofs()*I+s, NumLocalDofs()*J+s ) += mNa[I]*mNa[J]*MassFactor()*jac*qItr->Wt();
        }               
        
}

Numeric Formulation::SetDNx( list<Element>::const_iterator eItr, const Point &pt )
{
        nMatrix DNxi( eItr->NumNodes(), eItr->Dimension() );
        eItr->DNa( pt, DNxi );
        return grad_shape( mCoord, DNxi, mDNx );
}

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

} // end namespace

#endif

---