mapping.hh

#ifndef DUNE_FEM_MAPPING_HH
#define DUNE_FEM_MAPPING_HH
 
#include <iostream>
#include <vector>
 
namespace Dune
{
 
  namespace Fem
  {
 
 
    template<typename DFieldType,typename RFieldType, class DType, class RType>
    class Mapping
    {
    protected:
      typedef Mapping<DFieldType,RFieldType,DType,RType> MappingType;
 
    public:
      typedef DType DomainType;
 
      typedef RType  RangeType;
 
      typedef DFieldType DomainFieldType;
 
      typedef RFieldType RangeFieldType;
 
      Mapping()
        : lincomb_(),
          loopDetected_( false )
      {
        // store pointer to myself
        lincomb_.push_back( term( *this, 1. ) );
      }
 
      virtual ~Mapping ()
      {
      }
 
      MappingType& operator = (const MappingType &mapping)
      {
        const MappingType &m = dynamic_cast<const MappingType& >( mapping );
 
        lincomb_.erase( lincomb_.begin(), lincomb_.end() );
        typedef typename std::vector<term>::const_iterator iterator;
 
        // reserve memory
        lincomb_.reserve( m.lincomb_.size() );
 
        iterator end = m.lincomb_.end();
        for (iterator it = m.lincomb_.begin(); it != end; it++ )
        {
          lincomb_.push_back( *it );
        }
        return *this;
      }
 
      void operator() (const DomainType &arg, RangeType &dest ) const
      {
        int count = 0;
        typedef typename std::vector<term>::const_iterator const_iterator;
        const_iterator end = lincomb_.end();
        for ( const_iterator it = lincomb_.begin(); it != end; ++it )
        {
          const term& current = (*it);
          if ( count == 0 )
          {
            // call apply of term, see below
            current.apply( arg, dest );
          }
          else
          {
            // call applyAdd of term, see below
            current.applyAdd( arg, dest );
          }
          ++count;
        }
 
        loopDetected_ = false ;
      }
 
    private:
      virtual void apply (const DomainType &arg, RangeType &dest) const
      {
        // do noting if apply method leads to myself
        if( loopDetected_ )
        {
#ifndef NDEBUG
          std::cerr << "WARNING: Mapping::apply: Loop detected! Check overloading of operator() and apply! " << std::endl;
#endif
          return ;
        }
 
        // set applyLoop to true again
        loopDetected_ = true ;
 
        operator()(arg, dest);
      }
 
      struct term {
        term() : v_(NULL), scalar_(1.0), scaleIt_(false) { }
        term(const term& other)
          : v_(other.v_), scalar_(other.scalar_), scaleIt_(other.scaleIt_) { }
 
        term(const MappingType &mapping, RangeFieldType scalar ) : v_(&mapping), scalar_(scalar), scaleIt_( true ) {
          if ( scalar_ == 1. ) {
            scaleIt_ = false;
          }
        }
 
        void apply(const DomainType &arg, RangeType &dest) const
        {
          v_->apply( arg, dest );
          if ( scaleIt_ )
          {
            dest *= scalar_;
          }
        }
 
        void applyAdd(const DomainType &arg, RangeType &dest) const
        {
          // note, copying here might be costly
          RangeType tmp( dest );
 
          v_->apply( arg, tmp );
          if ( scalar_ == 1. )
          {
            dest += tmp;
          }
          else if ( scalar_ == -1. )
          {
            dest -= tmp;
          }
          else
          {
            tmp *= scalar_;
            dest += tmp;
          }
        }
 
      protected:
        const MappingType *v_;
        RangeFieldType scalar_;
        bool scaleIt_;
 
        // friendship for operations
        friend struct MappingOperators;
      };
 
 
      std::vector<term> lincomb_;
 
      mutable bool loopDetected_ ;
 
      // friendship for operations
      friend struct MappingOperators;
    };
 
 
    struct MappingOperators
    {
      template<typename DFieldType,typename RFieldType, class DType, class RType>
      static inline void
      copyMapping(const Mapping<DFieldType,RFieldType,DType,RType>& org,
                  Mapping<DFieldType,RFieldType,DType,RType>& copy)
      {
        typedef Mapping<DFieldType,RFieldType,DType,RType> MappingType;
        typedef typename std::vector< typename MappingType :: term > :: const_iterator iterator;
 
        // clear mapping entries
        copy.lincomb_.clear();
 
        {
          iterator end = org.lincomb_.end();
          for ( iterator it = org.lincomb_.begin(); it != end; ++it )
          {
            copy.lincomb_.push_back( *it );
          }
        }
      }
 
      template<typename DFieldType,typename RFieldType, class DType, class RType>
      static inline Mapping<DFieldType,RFieldType,DType,RType>
      addMappings(const Mapping<DFieldType,RFieldType,DType,RType>& a,
                  const Mapping<DFieldType,RFieldType,DType,RType>& b)
      {
        typedef Mapping<DFieldType,RFieldType,DType,RType> MappingType;
        // new mapping
        MappingType newMapping;
 
        // copy mapping
        copyMapping(a, newMapping);
 
        typedef typename std::vector< typename MappingType :: term > :: const_iterator iterator;
 
        iterator end = b.lincomb_.end();
        for ( iterator it = b.lincomb_.begin(); it != end; ++it )
        {
          newMapping.lincomb_.push_back( *it );
        }
        return newMapping;
      }
 
      template<typename DFieldType,typename RFieldType, class DType, class RType>
      static inline Mapping<DFieldType,RFieldType,DType,RType>
      substractMappings(const Mapping<DFieldType,RFieldType,DType,RType>& a,
                        const Mapping<DFieldType,RFieldType,DType,RType>& b)
      {
        typedef Mapping<DFieldType,RFieldType,DType,RType> MappingType;
        typedef typename  MappingType :: term term;
        typedef typename std::vector< term > :: const_iterator iterator;
        // new mapping
        MappingType newMapping;
 
        // copy mapping
        copyMapping(a, newMapping);
 
        iterator end = b.lincomb_.end();
        for ( iterator it = b.lincomb_.begin(); it != end; ++it )
        {
          newMapping.lincomb_.push_back( term( *it->v_, -it->scalar_ ) );
        }
        return newMapping;
      }
 
      template<typename DFieldType,typename RFieldType, class DType, class RType>
      static inline Mapping<DFieldType,RFieldType,DType,RType>
      multiplyMapping(const Mapping<DFieldType,RFieldType,DType,RType>& a,
                      const RFieldType& scalar)
      {
        typedef Mapping<DFieldType,RFieldType,DType,RType> MappingType;
        typedef typename  MappingType :: term term;
        typedef typename std::vector< term > :: iterator iterator;
        // new mapping
        MappingType newMapping;
 
        // copy mapping
        copyMapping(a, newMapping);
 
        iterator end = newMapping.lincomb_.end();
        for ( iterator it = newMapping.lincomb_.begin(); it != end; ++it )
        {
          it->scalar_ *= scalar;
        }
        return newMapping;
      }
 
      template<typename DFieldType,typename RFieldType, class DType, class RType>
      static inline Mapping<DFieldType,RFieldType,DType,RType>
      divideMapping(const Mapping<DFieldType,RFieldType,DType,RType>& a,
                    const RFieldType& scalar)
      {
        RFieldType factor = RFieldType(1)/scalar;
        return multiplyMapping(a,factor);
      }
    };
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator +(const Mapping<DFieldType,RFieldType,DType,RType>& a,
               const Mapping<DFieldType,RFieldType,DType,RType>& b)
    {
      return MappingOperators::addMappings(a,b);
    }
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator -(const Mapping<DFieldType,RFieldType,DType,RType>& a,
               const Mapping<DFieldType,RFieldType,DType,RType>& b)
    {
      return MappingOperators::substractMappings(a,b);
    }
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator *(const Mapping<DFieldType,RFieldType,DType,RType>& mapping,
               const RFieldType& factor)
    {
      return MappingOperators::multiplyMapping(mapping,factor);
    }
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator *(const RFieldType& factor,
               const Mapping<DFieldType,RFieldType,DType,RType>& mapping)
    {
      return MappingOperators::multiplyMapping(mapping,factor);
    }
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator /(const Mapping<DFieldType,RFieldType,DType,RType>& mapping,
               const RFieldType& factor)
    {
      return MappingOperators::divideMapping(mapping,factor);
    }
 
    template<class DFieldType, class RFieldType, class DType, class RType>
    static inline Mapping<DFieldType,RFieldType,DType,RType>
    operator /(const RFieldType& factor,
               const Mapping<DFieldType,RFieldType,DType,RType>& mapping)
    {
      return MappingOperators::divideMapping(mapping,factor);
    }
 
  } // namespace Fem
 
} // namespace Dune
#endif // #ifndef DUNE_FEM_MAPPING_HH