unityTriangleIntegratorShellFE.h

Go to the documentation of this file.

#ifndef __UNITTRIANGLEINTEGRATORSHELLFE_H
#define __UNITTRIANGLEINTEGRATORSHELLFE_H

namespace shellFE {

//=============================================================================
// Vector-valued Intefaces
//=============================================================================

template <typename ConfiguratorType, typename Imp>
class UnitTriangleFENonlinOpIntegratorShellFE  {

protected:
    typedef typename ConfiguratorType::RealType RealType;
    typedef typename ConfiguratorType::VectorType VectorType;
    typedef typename ConfiguratorType::ElementType ElementType;
    const ConfiguratorType & _config;

public:
    UnitTriangleFENonlinOpIntegratorShellFE( const ConfiguratorType & Conf ) : _config( Conf){ }

    virtual ~UnitTriangleFENonlinOpIntegratorShellFE( ) {}

    void assembleAdd ( VectorType &Dest ) const {
        RealType *nl_cache = new RealType[ _config.maxNumQuadPoints() ];

        for ( int elementIdx = 0; elementIdx < _config.getInitializer().getNumTriangs(); ++elementIdx){
            const ElementType& El ( _config.getInitializer().getTriang( elementIdx ) );
            const int numLocalDofs = _config.getNumLocalDofs ( El );

            const typename ConfiguratorType::BaseFuncSetType &bfs = _config.getBaseFunctionSet ( El );
            const int numQuadPoints = bfs.numQuadPoints( );

            for ( int q = 0; q < numQuadPoints; ++q )
                nl_cache[q] = this->asImp().getNonlinearity ( El, q );

            RealType aux;

            for ( int dof = 0; dof < numLocalDofs; ++dof ) {
                aux = 0.;

                for ( int q = 0; q < numQuadPoints; ++q )
                    aux += nl_cache[q] * bfs.evaluate ( dof, q ) * bfs.getWeight ( q );

                Dest[ _config.localToGlobal ( El, dof ) ] += 0.5 * aux;
            }
        }
        delete[] nl_cache;
    }

    RealType getNonlinearity ( const typename ConfiguratorType::ElementType & El, int QuadPoint ) const {
        throw std::invalid_argument( aol::strprintf ( "Called the interface function. In File %s at line %d.", __FILE__, __LINE__ ).c_str() );
        return this->asImp().getNonlinearity ( El, QuadPoint );
    }

protected:
    // barton-nackman
    inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
    inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }

};

//=============================================================================
// Matrix-valued Intefaces
//=============================================================================


template < typename ConfiguratorType, typename Imp >
class MatrixValuedIntegratorBase {
public:
  typedef typename ConfiguratorType::RealType   RealType;
  typedef typename ConfiguratorType::LocalMatrixType LocalMatrixType;
  typedef typename ConfiguratorType::TripletType TripletType;
  typedef typename ConfiguratorType::ElementType ElementType;
  typedef typename ConfiguratorType::MaskType MaskType;

  explicit MatrixValuedIntegratorBase ( const ConfiguratorType &conf ) : _config ( conf ) {}

protected:
  void assembleTripletList ( std::vector<TripletType> & tripletList, const RealType Factor ) const {
    tripletList.reserve(_config.getInitializer().getNumTriangs() * aol::Sqr( _config.getNumLocalDofs() ) );
    LocalMatrixType localMatrix;
    int globalDofs[ ConfiguratorType::maxNumLocalDofs ];
    for ( int elementIdx = 0; elementIdx < _config.getInitializer().getNumTriangs(); ++elementIdx){
      const ElementType& El ( _config.getInitializer().getTriang( elementIdx ) );
      this->asImp().prepareLocalMatrix ( El, localMatrix );
      const int numLocalDofs = _config.getNumLocalDofs ( El );

      // Part of exercise 1
    }
  }
  
public:

  template <typename SparseMatrixType>
  void assemble ( SparseMatrixType &Dest, const RealType Factor = 1.0 ) const {
    std::vector<TripletType> tripletList;
    assembleTripletList ( tripletList, Factor );
    Dest.setFromTriplets( tripletList.cbegin(), tripletList.cend() );
  }
  
  template <typename SparseMatrixType>
  void assembleDirichlet ( SparseMatrixType &Dest, const MaskType& boundaryMask, const RealType Factor = 1.0 ) const {
       
    std::vector<TripletType> tripletList;
    assembleTripletList ( tripletList, Factor );
    
    std::vector<TripletType> tripletListMasked;
    tripletListMasked.reserve(_config.getInitializer().getNumTriangs() * aol::Sqr( _config.getNumLocalDofs() ) );

    for( unsigned iter=0; iter < tripletList.size(); ++iter ){
      if( (boundaryMask[tripletList[iter].row()]) || (boundaryMask[tripletList[iter].col()]) ){
       //Boundary node!        
      } else {
        tripletListMasked.push_back( tripletList[iter] );
      }
    }
    
    for ( int i = 0; i < _config.getNumGlobalDofs(); ++i ){
      if ( boundaryMask[i] )
         tripletListMasked.push_back( TripletType( i, i, 1.0 ) );
    }
    
    Dest.setFromTriplets( tripletListMasked.begin(), tripletListMasked.end() );
  }
  
protected:
  // barton-nackman
  inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
  inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
  
  const ConfiguratorType &_config;
};




template <typename ConfiguratorType, typename Imp >
class UnitTriangleFELinWeightedStiffIntegrator :
      public MatrixValuedIntegratorBase< ConfiguratorType, UnitTriangleFELinWeightedStiffIntegrator<ConfiguratorType, Imp> > {
protected:
  typedef typename ConfiguratorType::RealType RealType;
  typedef typename ConfiguratorType::DomVecType DomVecType;
  typedef typename ConfiguratorType::Matrix22 Matrix22;
  typedef typename ConfiguratorType::LocalMatrixType LocalMatrixType;
  const ConfiguratorType &_config;
  
public:
  UnitTriangleFELinWeightedStiffIntegrator ( const ConfiguratorType & Config ) :
   MatrixValuedIntegratorBase<  ConfiguratorType, UnitTriangleFELinWeightedStiffIntegrator<ConfiguratorType, Imp> > ( Config ),
  _config ( Config ) {}


  inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint,
                               Matrix22 &Matrix ) const {
    this->asImp().getCoeffMatrix ( El, QuadPoint, Matrix );
  }

  inline void prepareLocalMatrix ( const typename ConfiguratorType::ElementType &El, 
                                  LocalMatrixType &LocalMatrix ) const {
    const int numDofs = _config.getNumLocalDofs ( El ); 
    
    for ( int i = 0; i < numDofs; ++i )
      for ( int j = 0; j < numDofs; ++j )
        LocalMatrix(i,j) = 0.;

    Matrix22 mat;
    DomVecType matgrad1;

    const typename ConfiguratorType::BaseFuncSetType &bfs = _config.getBaseFunctionSet ( El );
    const int numQuadPoints = bfs.numQuadPoints( );

    // Part of exercise 1
  }

protected:
  inline Imp &asImp() { return static_cast<Imp&> ( *this ); }
  inline const Imp &asImp() const { return static_cast<const Imp&> ( *this ); }

};

} // end namespace

#endif