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

Go to the documentation of this file.

#ifndef _SOLVER_BASSO_H_
#define _SOLVER_BASSO_H_

#include "basso.h"

namespace basso {

template < class KMAT >
class Solver {

public:
        
        typedef KMAT matrix_type;

        Solver( KMAT &K, nArray &f )  
                { Kptr=&K; fptr=&f; 
                        wt=mat_trace(K)/static_cast<Numeric>(f.size()); 
                        solved=false; 
                }
        
        int NumEquations() const { return fptr->size(); }
        
        virtual int Solve( nArray &d ) = 0;
        
        template < class INDX >
        int SolveConstrained( nArray &d, INDX &ifix );
        
        int SolveConstrained( nArray &d, iArray &ifix, nArray &ival );

        virtual void GetReactionForces( const nArray &d, nArray &freac )
                { mult( Kreac, d, freac ); }

        bool Solved() const { return solved; }
        
protected:
        
        template < class INDX >
        void GetReactionMatrix( const INDX &ifix );
        
        template < class INDX >
        void SetConstraintsOnMatrix( const INDX &ifix );
        
        template < class INDX >
        void SetConstraintsOnF( const INDX &ifix );
        
        template < class INDX >
        void SetConstraintsOnF( const INDX &ifix, nArray &ival );

protected:
        KMAT *Kptr;
        nArray *fptr;
        KMAT Kreac;
        Numeric wt;
        bool solved;
};


template < class KMAT> template< class INDX >
void Solver<KMAT>::GetReactionMatrix( const INDX &ifix )  
{
        int nc=ifix.size(), ndof=mat_nrows(*Kptr);
         
        // get reaction equations
        if ( mat_nrows(Kreac)!=nc || mat_ncols(Kreac)!=ndof )
                resize( Kreac, nc, ndof );
        typename INDX::const_iterator ifixItr;
        int i=0;
        for ( ifixItr=ifix.begin(); ifixItr!=ifix.end(); ++ifixItr, ++i ) {
                int ii=*ifixItr;
                copy ( sub_matrix(*Kptr,sub_interval(ii,1),sub_interval(0,ndof)), 
                                   sub_matrix(Kreac,sub_interval(i,1),  sub_interval(0,ndof)) );
        }       
}

template < class KMAT> template< class INDX >
void Solver<KMAT>::SetConstraintsOnMatrix( const INDX &ifix )  
{

        typename INDX::const_iterator ifixItr;
        for ( ifixItr=ifix.begin(); ifixItr!=ifix.end(); ++ifixItr )
        {
                int ii=*ifixItr;
                (*Kptr)(ii,ii)=wt;
                for ( int jj=0; jj<NumEquations(); ++jj )
                        if ( ii != jj)
                        {
                                (*Kptr)(ii,jj)=0.0; 
                                (*Kptr)(jj,ii)=0.0;
                        }
        }       
}

template < class KMAT> template< class INDX >
void Solver<KMAT>::SetConstraintsOnF( const INDX &ifix )  
{
        typename INDX::const_iterator ifixItr;
        for ( ifixItr=ifix.begin(); ifixItr!=ifix.end(); ++ifixItr )
                (*fptr)[*ifixItr]=0.0; 
}

template < class KMAT> template< class INDX >
void Solver<KMAT>::SetConstraintsOnF( const INDX &ifix, nArray &ival )  
{
        // modify  rhs f := -K(:,ifix)*ival + f 
        mult_add( sub_matrix( *Kptr, sub_interval(0,NumEquations()), sub_index(ifix) ), 
                scaled(ival,-1.0), *fptr );
        typename INDX::const_iterator ifixItr;
        int i=0;
        for ( ifixItr=ifix.begin(); ifixItr!=ifix.end(); ++ifixItr, ++i )
                (*fptr)[*ifixItr]=ival[i]*wt; 
}

template < class KMAT> template< class INDX >
int Solver<KMAT>::SolveConstrained( nArray &d, INDX &ifix ) 
{
        GetReactionMatrix( ifix );
        SetConstraintsOnMatrix( ifix );
        SetConstraintsOnF( ifix );
        return Solve(d); 
}

template < class KMAT>
int Solver<KMAT>::SolveConstrained( nArray &d, iArray &ifix, nArray &ival )
{
        GetReactionMatrix( ifix );
        SetConstraintsOnMatrix( ifix, ival );
        SetConstraintsOnF( ifix, ival );
        return Solve(d);        
}

} // end namespace

#endif

---