localmatrix.hh

#ifndef DUNE_FEM_LOCALMATRIX_HH
#define DUNE_FEM_LOCALMATRIX_HH
 
#include <optional>
 
//- Dune includes
#include <dune/fem/misc/bartonnackmaninterface.hh>
#include "../../common/explicitfieldvector.hh"
 
namespace Dune
{
 
  namespace Fem
  {
 
    //- forward declaration of MatrixColumnObject, see below
    template <class Traits>
    class MatrixColumnObject ;
 
 
 
    // LocalMatrixInterface
    // --------------------
 
    template< class LocalMatrixTraits >
    class LocalMatrixInterface
    : public BartonNackmanInterface< LocalMatrixInterface< LocalMatrixTraits >,
                                     typename LocalMatrixTraits::LocalMatrixType >
    {
      typedef LocalMatrixInterface< LocalMatrixTraits > ThisType;
      typedef BartonNackmanInterface< LocalMatrixInterface< LocalMatrixTraits >,
                                      typename LocalMatrixTraits::LocalMatrixType >
        BaseType;
 
    public:
      typedef LocalMatrixTraits Traits;
 
      typedef ThisType LocalMatrixInterfaceType;
 
      typedef typename Traits :: LocalMatrixType LocalMatrixType;
 
      typedef typename Traits :: RangeFieldType RangeFieldType;
 
      typedef typename Traits :: DomainSpaceType DomainSpaceType;
 
      typedef typename Traits :: RangeSpaceType RangeSpaceType;
 
      typedef typename DomainSpaceType :: BasisFunctionSetType
        DomainBasisFunctionSetType;
 
      typedef typename RangeSpaceType :: BasisFunctionSetType
        RangeBasisFunctionSetType;
 
      typedef typename DomainSpaceType::EntityType DomainEntityType;
      typedef typename RangeSpaceType::EntityType RangeEntityType;
 
      typedef typename Traits :: LittleBlockType  LittleBlockType;
 
      typedef MatrixColumnObject< Traits >  MatrixColumnType;
 
    protected:
      using BaseType::asImp;
 
      LocalMatrixInterface ()
      {}
 
    public:
      void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION
          ( asImp().init( domainEntity, rangeEntity ) );
      }
 
      void bind ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION
          ( asImp().bind( domainEntity, rangeEntity ) );
      }
 
      void unbind ()
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().unbind() );
      }
 
      void add ( const int localRow,
                 const int localCol,
                 const RangeFieldType &value )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
          asImp().add(localRow,localCol,value));
      }
 
      void set ( const int localRow,
                 const int localCol,
                 const RangeFieldType &value )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
          asImp().set(localRow,localCol,value));
      }
 
      void clearRow( const int  localRow )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
            asImp().clearRow( localRow ));
      }
 
      void clearCol( const int localCol )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
            asImp().clearCol( localCol ));
      }
 
 
      template <class DomainLocalFunctionType,
                class RangeLocalFunctionType>
      void multiplyAdd(const DomainLocalFunctionType& lhs,
                       RangeLocalFunctionType& rhs) const
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
          asImp().multiplyAdd( lhs, rhs ) );
      }
 
      const RangeFieldType get ( const int localRow,
                                 const int localCol ) const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().get(localRow,localCol));
        return asImp().get(localRow,localCol);
      }
 
      void scale ( const RangeFieldType& scalar )
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(
          asImp().scale( scalar ) );
      }
 
      void clear ()
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().clear());
      }
 
      void resort ()
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().resort());
      }
 
      int rows () const
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().rows() );
        return asImp().rows();
      }
 
      int columns () const
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().columns() );
        return asImp().columns();
      }
 
      const DomainSpaceType &domainSpace () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().domainSpace() );
        return asImp().domainSpace();
      }
 
      const RangeSpaceType &rangeSpace () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().rangeSpace() );
        return asImp().rangeSpace();
      }
 
      const DomainBasisFunctionSetType &domainBasisFunctionSet () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().domainBasisFunctionSet() );
        return asImp().domainBasisFunctionSet();
      }
 
      const RangeBasisFunctionSetType &rangeBasisFunctionSet () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().rangeBasisFunctionSet() );
        return asImp().rangeBasisFunctionSet();
      }
 
      const DomainEntityType &domainEntity () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().domainEntity() );
        return asImp().domainEntity();
      }
 
      const RangeEntityType &rangeEntity () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( asImp().rangeEntity() );
        return asImp().rangeEntity();
      }
 
      MatrixColumnType column( const unsigned int col )
      {
        return MatrixColumnType( asImp(), col );
      }
 
      void finalize()
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().finalize());
      }
    };
 
 
 
    // LocalMatrixDefault
    // ------------------
 
    template< class LocalMatrixTraits >
    class LocalMatrixDefault
    : public LocalMatrixInterface< LocalMatrixTraits >
    {
      typedef LocalMatrixDefault< LocalMatrixTraits > ThisType;
      typedef LocalMatrixInterface< LocalMatrixTraits > BaseType;
 
    public:
      typedef LocalMatrixTraits Traits;
 
      typedef typename BaseType::DomainSpaceType DomainSpaceType;
      typedef typename BaseType::RangeSpaceType RangeSpaceType;
 
      typedef typename BaseType::DomainBasisFunctionSetType DomainBasisFunctionSetType;
      typedef typename BaseType::RangeBasisFunctionSetType RangeBasisFunctionSetType;
 
      typedef typename BaseType::DomainEntityType DomainEntityType;
      typedef typename BaseType::RangeEntityType RangeEntityType;
 
    protected:
      const DomainSpaceType &domainSpace_;
      const RangeSpaceType &rangeSpace_;
 
      DomainBasisFunctionSetType domainBaseSet_;
      RangeBasisFunctionSetType rangeBaseSet_;
 
      std::optional< DomainEntityType > domainEntity_;
      std::optional< RangeEntityType > rangeEntity_;
 
    protected:
      LocalMatrixDefault ( const DomainSpaceType &domainSpace,
                           const RangeSpaceType &rangeSpace )
      : domainSpace_( domainSpace ),
        rangeSpace_( rangeSpace ),
        domainBaseSet_(),
        rangeBaseSet_()
      {}
 
      template< class DomainEntityType, class RangeEntityType >
      LocalMatrixDefault ( const DomainSpaceType &domainSpace,
                           const RangeSpaceType &rangeSpace,
                           const DomainEntityType &domainEntity,
                           const RangeEntityType &rangeEntity )
      : domainSpace_( domainSpace ),
        rangeSpace_( rangeSpace ),
        domainBaseSet_(),
        rangeBaseSet_()
      {
        bind( domainEntity, rangeEntity );
      }
 
      LocalMatrixDefault ( const LocalMatrixDefault& org )
      : domainSpace_( org.domainSpace_ ),
        rangeSpace_( org.rangeSpace_ ),
        domainBaseSet_( org.domainBaseSet_ ),
        rangeBaseSet_( org.rangeBaseSet_ ),
        domainEntity_( org.domainEntity_ ),
        rangeEntity_( org.rangeEntity_ )
      {}
 
    public:
      void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        bind( domainEntity, rangeEntity );
      }
 
      void bind ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        domainBaseSet_ = domainSpace_.basisFunctionSet( domainEntity );
        rangeBaseSet_ = rangeSpace_.basisFunctionSet( rangeEntity );
        domainEntity_.emplace(domainEntity);
        rangeEntity_.emplace(rangeEntity);
      }
 
      void unbind ()
      {
        domainBaseSet_ = DomainBasisFunctionSetType();
        rangeBaseSet_  = RangeBasisFunctionSetType();
        domainEntity_.reset();
        rangeEntity_.reset();
      }
 
      void resort () {}
 
      void finalize () {}
 
      int rows () const { return rangeBaseSet_.size(); }
 
      int columns () const { return domainBaseSet_.size(); }
 
      const DomainSpaceType &domainSpace () const { return domainSpace_; }
 
      const RangeSpaceType &rangeSpace () const { return rangeSpace_; }
 
      const DomainBasisFunctionSetType &domainBasisFunctionSet () const
      {
        return domainBaseSet_;
      }
 
      const RangeBasisFunctionSetType &rangeBasisFunctionSet () const
      {
        return rangeBaseSet_;
      }
 
      const DomainEntityType &domainEntity () const { return domainEntity_.value(); }
      const RangeEntityType &rangeEntity () const { return rangeEntity_.value(); }
 
      template <class DomainLocalFunctionType,
                class RangeLocalFunctionType>
      void multiplyAdd(const DomainLocalFunctionType& lhs,
                              RangeLocalFunctionType& rhs) const
      {
        const int row = this->rows();
        const int col = this->columns();
        for(int i=0; i<row; ++i)
        {
          for(int j=0; j<col; ++j)
          {
            rhs[i] += this->get(i,j) * lhs[j];
          }
        }
      }
 
      void clearRow( const int localRow )
      {
        const int col = this->columns();
        for(int j = 0; j < col; ++j)
        {
          this->set(localRow, j, 0);
        }
      }
 
      void clearCol( const int localCol )
      {
        const int row = this->rows();
        for(int i = 0; i < row; ++i)
        {
          this->set(i, localCol, 0);
        }
      }
    };
 
    template <class Traits>
    class MatrixColumnObject
    {
    public:
      typedef typename Traits :: LocalMatrixType LocalMatrixType;
 
      typedef typename Traits :: RangeSpaceType  RangeSpaceType;
 
      typedef typename RangeSpaceType :: RangeType          RangeType ;
      typedef typename RangeSpaceType :: JacobianRangeType  JacobianRangeType ;
      typedef typename RangeSpaceType :: RangeFieldType     RangeFieldType ;
 
    protected:
      // reference to local matrix
      LocalMatrixType& localMatrix_;
      // local column number
      const unsigned int column_;
 
      MatrixColumnObject( LocalMatrixType& localMatrix, const unsigned int col )
        : localMatrix_( localMatrix ),
          column_( col )
      {
      }
 
      // at the moment only allow LocalMatrixInterface to construct this object
      friend class LocalMatrixInterface< Traits >;
 
    public:
 
      template <class RangeVectorType>
      void axpy( const RangeVectorType& phi,
                 const Explicit<RangeType>& factor,
                 const RangeFieldType& weight = RangeFieldType(1) )
      {
        const unsigned int numBasisFunctions = localMatrix_.rows();
        assert( phi.size() >= numBasisFunctions );
        for( unsigned int row = 0; row < numBasisFunctions; ++ row )
        {
          RangeFieldType value = factor * phi[ row ];
          localMatrix_.add( row, column_, weight * value );
        }
      }
 
      template <class JacobianVectorType>
      void axpy( const JacobianVectorType& dphi,
                 const JacobianRangeType& jacobianFactor,
                 const RangeFieldType& weight = RangeFieldType(1) )
      {
        const unsigned int numBasisFunctions = localMatrix_.rows();
        assert( dphi.size() >= numBasisFunctions );
        for( unsigned int row = 0; row < numBasisFunctions; ++ row )
        {
          RangeFieldType value = 0;
          for( int k = 0; k < jacobianFactor.rows; ++k )
            value += jacobianFactor[ k ] * dphi[ row ][ k ];
 
          localMatrix_.add( row, column_, weight * value );
        }
      }
 
      template <class RangeVectorType, class JacobianVectorType>
      void axpy( const RangeVectorType& phi,
                 const JacobianVectorType& dphi,
                 const Explicit<RangeType>& factor,
                 const JacobianRangeType& jacobianFactor,
                 const RangeFieldType& weight = RangeFieldType(1) )
      {
        const unsigned int numBasisFunctions = localMatrix_.rows();
        assert( phi.size() >= numBasisFunctions );
        assert( dphi.size() >= numBasisFunctions );
        for( unsigned int row = 0; row < numBasisFunctions; ++ row )
        {
          RangeFieldType value = factor * phi[ row ];
          for( int k = 0; k < jacobianFactor.rows; ++k )
            value += jacobianFactor[ k ] * dphi[ row ][ k ];
 
          localMatrix_.add( row, column_, weight * value );
        }
      }
    };
 
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_LOCALMATRIX_HH