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

Go to the documentation of this file.

#ifndef _TRACTION_BC_BASSO_H_
#define _TRACTION_BC_BASSO_H_

// std includes
#include <list>

// Basso includes
#include "basso.h"
#include "NeumannBC.h"
#include "basic_fem_operations.h"
#include "geometry.h"

namespace basso {


class TractionBC : public NeumannBC {

public:
        
        // PUBLIC DEFS
        
        // CONSTRUCTORS 
        TractionBC ( list< Element > &elem, const Array< Node > &nodes, 
                const Array<Dof> &dofs, const nArray &n, Numeric sf=1.0 ) : NeumannBC ( dofs, sf )
                { fElement=&elem; fNode=&nodes; mDirection=n; }

        TractionBC ( list< Element > &elem, const Array< Node > &nodes, 
                const Dof &ldof, Numeric sf=1.0 ) : NeumannBC ( ldof, sf ), mDirection(1)
                { fElement=&elem; fNode=&nodes; mDirection[0]=1.0; }

        
        virtual ~TractionBC() {  }
        
        // ACCESSORS
        
        // MEMBER FUNCTIONS
        virtual void ContributeForce( const DofMap &dofmap, nArray &f );
        
protected:
        virtual Numeric Jacobian( list<Element>::const_iterator eItr, const Point &pt );
        virtual void GetElementForce( list<Element>::const_iterator eItr, nArray &fe );
        virtual void GetElementForceLine2( list<Element>::const_iterator eItr, nArray &fe );
        virtual void GetElementForceTria3( list<Element>::const_iterator eItr, nArray &fe );
        virtual bool SetElementParameters( list<Element>::const_iterator eItr, int qorder );

        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(); }

protected:
        
        nArray mDirection;
        list< Element > *fElement;
        const Array< Node > *fNode;

        nMatrix mCoord, mDNxi;
        nArray mNa;
        int mNne, mSdim, mEdim;
        QuadratureRule mQuadrature;
        
        BasisType mLastElemType;
};


void TractionBC::GetElementForceLine2( list<Element>::const_iterator eItr, nArray &fe )
{
        int n0=eItr->NodeId(0), n1=eItr->NodeId(1);
        Numeric le = dist( Point( (*(this->fNode))[n1] ), Point( (*(this->fNode))[n0] ) );
        int i=0;
        for ( int n=0; n<2; ++n )
                for ( int d=0; d<this->mDofs.size(); ++d )
                        fe[i++] = 0.5*le*(this->mScaleFactor)*(this->mDirection[d]);
}

void TractionBC::GetElementForceTria3( list<Element>::const_iterator eItr, nArray &fe )
{
        int n0=eItr->NodeId(0), n1=eItr->NodeId(1), n2=eItr->NodeId(2);
        Numeric ae = area_triangle( Point( (*(this->fNode))[n0] ), 
                Point( (*(this->fNode))[n1] ), Point( (*(this->fNode))[n2] ) );
        int i=0;
        for ( int n=0; n<3; ++n )
                for ( int d=0; d<this->mDofs.size(); ++d )
                        fe[i++] = 0.3333333333333333333333*ae*(this->mScaleFactor)*(this->mDirection[d]);
}

void TractionBC::GetElementForce( list<Element>::const_iterator eItr, nArray &fe )
{

        int fdim = eItr->NumNodes()*(this->mDofs.size());
        resize( fe, fdim );
        clear(fe);
        
        if ( eItr->Type() == kLINE2 )
        {
                GetElementForceLine2( eItr, fe );
                return;
        }
        else if ( eItr->Type() == kTRIA3 )
        {
                GetElementForceTria3( eItr, fe );
                return;
        }
        
        SetElementParameters( eItr, eItr->Order() );    
        Numeric jac;
        QuadratureRule::const_iterator qItr;
        for ( qItr=QuadratureBegin(); qItr!=QuadratureEnd(); ++qItr )   // quadrature loop
        {
                jac = Jacobian( eItr, qItr->Pt() );
                int i=0;
                for ( int n=0; n<mNne; ++n )
                        for ( int d=0; d<this->mDofs.size(); ++d )
                                fe[i++] += (this->mNa[n])*(this->mScaleFactor)*(this->mDirection[d]);
        }
        
}

void TractionBC::ContributeForce( const DofMap &dofmap, nArray &f )
{
        nArray fe;
        list< Element >::const_iterator eItr;
        Element::const_iterator nItr;
 
        for ( eItr=ElementBegin(); eItr!=ElementEnd(); ++eItr )  // element loop
        {
                GetElementForce( eItr, fe );
                dofmap.Scatter( eItr->Connectivity(), this->mDofs, fe, f );
        }
}

Numeric TractionBC::Jacobian( list<Element>::const_iterator eItr, const Point &pt )
{ 
        eItr->DNa( pt, mDNxi );
        if ( mEdim==1 )
                return jacobian13( mCoord, mDNxi );
        else            
                return jacobian23( mCoord, mDNxi );
}

bool TractionBC::SetElementParameters( list<Element>::const_iterator eItr, int qorder )
{
        mSdim = 3;
        mEdim = 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, mEdim ); 
        resize( mNa, mNne );
        resize( mCoord, mNne, mSdim ); 
        eItr->Quadrature( qorder, mQuadrature );
        mLastElemType = eItr->Type(); 
        return true;
}


} // end namespace

#endif

---