assemblerutilities.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
#define DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
 
#include <algorithm>
#include <tuple>
 
#include <dune/pdelab/constraints/common/constraintstransformation.hh>
#include <dune/pdelab/gridoperator/common/localmatrix.hh>
 
namespace Dune{
  namespace PDELab{
 
    template<typename GO>
    struct LocalAssemblerTraits
    {
 
      typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
 
      typedef typename GO::Traits::TestGridFunctionSpace TestGridFunctionSpace;
 
 
      typedef typename GO::Traits::TrialGridFunctionSpaceConstraints TrialGridFunctionSpaceConstraints;
 
      typedef typename GO::Traits::TestGridFunctionSpaceConstraints TestGridFunctionSpaceConstraints;
 
 
      typedef typename GO::Traits::MatrixBackend MatrixBackend;
 
 
      typedef typename GO::Traits::DomainField DomainField;
 
      typedef typename GO::Traits::Domain Domain;
 
      typedef typename GO::Traits::Domain Solution;
 
 
      typedef typename GO::Traits::RangeField RangeField;
 
      typedef typename GO::Traits::Range Range;
 
      typedef typename GO::Traits::Range Residual;
 
 
      typedef typename GO::Traits::JacobianField JacobianField;
 
      typedef typename GO::Traits::Jacobian Jacobian;
 
      typedef typename MatrixBackend::template Pattern<
        Jacobian,
        TestGridFunctionSpace,
        TrialGridFunctionSpace
        > MatrixPattern;
 
      typedef typename GO::BorderDOFExchanger BorderDOFExchanger;
 
    };
 
    // ********************************************************************************
    // default local pattern implementation
    // ********************************************************************************
 
    class SparsityLink : public std::tuple<int,int>
    {
    public:
      SparsityLink ()
      {}
 
      SparsityLink (int i, int j)
        : std::tuple<int,int>(i,j)
      {}
 
      inline int i () const
      {
        return std::get<0>(*this);
      }
 
      inline int j () const
      {
        return std::get<1>(*this);
      }
 
      void set (int i, int j)
      {
        std::get<0>(*this) = i;
        std::get<1>(*this) = j;
      }
    };
 
    class LocalSparsityPattern
      : public std::vector<SparsityLink>
    {
 
      // make push_back() inaccessible
      using std::vector<SparsityLink>::push_back;
 
    public:
 
 
      template<typename LFSV, typename LFSU>
      void addLink(const LFSV& lfsv, std::size_t i, const LFSU& lfsu, std::size_t j)
      {
        std::vector<SparsityLink>::push_back(
          SparsityLink(
            lfsv.localIndex(i),
            lfsu.localIndex(j)
          )
        );
      }
 
    };
 
 
 
    // ********************************************************************************
    // Assembler base class
    // ********************************************************************************
 
    template<typename B,
             typename CU=EmptyTransformation,
             typename CV=EmptyTransformation>
    class LocalAssemblerBase
    {
    public:
 
      typedef typename B::size_type SizeType;
 
      LocalAssemblerBase ()
        : pconstraintsu(&emptyconstraintsu), pconstraintsv(&emptyconstraintsv)
      {}
 
      LocalAssemblerBase (const CU& cu, const CV& cv)
        :pconstraintsu(&cu), pconstraintsv(&cv)
      {}
 
      const CU& trialConstraints() const
      {
        return *pconstraintsu;
      }
 
      const CV& testConstraints() const
      {
        return *pconstraintsv;
      }
 
 
      template<typename X>
      typename std::enable_if<
        AlwaysTrue<X>::value && !std::is_same<
          CV,
          EmptyTransformation
          >::value
        >::type
      forwardtransform(X & x, const bool postrestrict = false) const
      {
        typedef typename CV::const_iterator global_col_iterator;
        for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
          typedef typename global_col_iterator::value_type::first_type GlobalIndex;
          const GlobalIndex & contributor = cit->first;
 
          typedef typename global_col_iterator::value_type::second_type ContributedMap;
          typedef typename ContributedMap::const_iterator global_row_iterator;
          const ContributedMap & contributed = cit->second;
          global_row_iterator it  = contributed.begin();
          global_row_iterator eit = contributed.end();
 
          for(;it!=eit;++it)
            {
              // typename X::block_type block(x[contributor]);
              // block *= it->second;
              // x[it->first] += block;
              x[it->first] += it->second * x[contributor];
            }
        }
 
        if(postrestrict)
          for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit)
            x[cit->first]=0.;
      }
 
 
      // Disable forwardtransform for EmptyTransformation
      template<typename X>
      typename std::enable_if<
        AlwaysTrue<X>::value && std::is_same<
          CV,
          EmptyTransformation
          >::value
        >::type
      forwardtransform(X & x, const bool postrestrict = false) const
      {}
 
 
      template<typename X>
      typename std::enable_if<
        AlwaysTrue<X>::value && !std::is_same<
          CV,
          EmptyTransformation
          >::value
        >::type
      backtransform(X & x, const bool prerestrict = false) const
      {
        typedef typename CV::const_iterator global_col_iterator;
        for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
          typedef typename global_col_iterator::value_type::first_type GlobalIndex;
          const GlobalIndex & contributor = cit->first;
 
          typedef typename global_col_iterator::value_type::second_type ContributedMap;
          typedef typename ContributedMap::const_iterator global_row_iterator;
          const ContributedMap & contributed = cit->second;
          global_row_iterator it  = contributed.begin();
          global_row_iterator eit = contributed.end();
 
          if(prerestrict)
            x[contributor] = 0.;
 
          for(;it!=eit;++it)
            {
              // typename X::block_type block(x[it->first]);
              // block *= it->second;
              // x[contributor] += block;
              x[contributor] += it->second * x[it->first]; // PB: 27 Sep 12 this was the old version
            }
        }
      }
 
      // disable backtransform for empty transformation
      template<typename X>
      typename std::enable_if<
        AlwaysTrue<X>::value && std::is_same<
          CV,
          EmptyTransformation
          >::value
        >::type
      backtransform(X & x, const bool prerestrict = false) const
      {}
 
 
    protected:
 
      template<typename GCView, typename T>
      void eread (const GCView& globalcontainer_view, LocalMatrix<T>& localcontainer) const
      {
        for (int i = 0; i < localcontainer.N(); ++i)
          for (int j = 0; j < localcontainer.M(); ++j)
            localcontainer(i,j) = globalcontainer_view(i,j);
      }
 
      template<typename T, typename GCView>
      void ewrite (const LocalMatrix<T>& localcontainer, GCView& globalcontainer_view) const
      {
        for (int i = 0; i < localcontainer.N(); ++i)
          for (int j = 0; j < localcontainer.M(); ++j)
            globalcontainer_view(i,j) = localcontainer(i,j);
      }
 
      template<typename T, typename GCView>
      void eadd (const LocalMatrix<T>& localcontainer, GCView& globalcontainer_view) const
      {
        for (int i = 0; i < localcontainer.N(); ++i)
          for (int j = 0; j < localcontainer.M(); ++j)
            globalcontainer_view.add(i,j,localcontainer(i,j));
      }
 
      template<typename M, typename GCView>
      typename std::enable_if<
        AlwaysTrue<M>::value && !std::is_same<
          CV,
          EmptyTransformation
          >::value
      >::type
      scatter_jacobian(M& local_container, GCView& global_container_view, bool symmetric_mode) const
      {
        typedef typename GCView::RowIndexCache LFSVIndexCache;
        typedef typename GCView::ColIndexCache LFSUIndexCache;
 
        const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
        const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
 
        if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
          if (symmetric_mode)
            etadd_symmetric(local_container,global_container_view);
          else
            etadd(local_container,global_container_view);
        else
          {
 
            typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
            const LFSV& lfsv = lfsv_indices.localFunctionSpace();
 
            typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
            const LFSU& lfsu = lfsu_indices.localFunctionSpace();
 
            // optionally clear out columns that belong to Dirichlet-constrained DOFs to keep matrix symmetric
            if (symmetric_mode)
              {
                typedef typename LFSUIndexCache::ContainerIndex CI;
 
                for (size_t j = 0; j < lfsu_indices.size(); ++j)
                  {
                    const CI& container_index = lfsu_indices.containerIndex(j);
                    const typename CU::const_iterator cit = pconstraintsu->find(container_index);
                    if (cit != pconstraintsu->end())
                      {
                        // make sure we only have Dirichlet constraints
                        assert(cit->second.empty());
                        // clear out the current column
                        for (size_t i = 0; i < lfsv_indices.size(); ++i)
                          {
                            // we do not need to update the residual, since the solution
                            // (i.e. the correction) for the Dirichlet DOF is 0 by definition
                            local_container(lfsv,i,lfsu,j) = 0.0;
                          }
                      }
                  }
              }
 
            // write entries without considering constraints.
            // Dirichlet-constrained rows will be fixed in a postprocessing step.
            for (auto it = local_container.begin(); it != local_container.end(); ++it)
              {
                // skip 0 entries because they might not be present in the pattern
                if (*it == 0.0)
                  continue;
                global_container_view.add(it.row(),it.col(),*it);
              }
          }
      }
 
      // specialization for empty constraints container
      template<typename M, typename GCView>
      typename std::enable_if<
        AlwaysTrue<M>::value && std::is_same<
          CV,
          EmptyTransformation
          >::value
      >::type
      scatter_jacobian(M& local_container, GCView& global_container_view, bool symmetric_mode) const
      {
        // write entries without considering constraints.
        // Dirichlet-constrained rows will be fixed in a postprocessing step.
        for (auto it = local_container.begin(); it != local_container.end(); ++it)
          {
            // skip 0 entries because they might not be present in the pattern
            if (*it == 0.0)
              continue;
            global_container_view.add(it.row(),it.col(),*it);
          }
      }
 
      template<typename M, typename GCView>
      void etadd_symmetric (M& localcontainer, GCView& globalcontainer_view) const
      {
 
        typedef typename GCView::RowIndexCache LFSVIndexCache;
        typedef typename GCView::ColIndexCache LFSUIndexCache;
 
        const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
        const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
 
        typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
        const LFSV& lfsv = lfsv_indices.localFunctionSpace();
 
        typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
        const LFSU& lfsu = lfsu_indices.localFunctionSpace();
 
        for (size_t j = 0; j < lfsu_indices.size(); ++j)
          {
            if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
              {
                // clear out the current column
                for (size_t i = 0; i < lfsv_indices.size(); ++i)
                  {
                    // we do not need to update the residual, since the solution
                    // (i.e. the correction) for the Dirichlet DOF is 0 by definition
                    localcontainer(lfsv,i,lfsu,j) = 0.0;
                  }
              }
          }
 
        // hand off to standard etadd() method
        etadd(localcontainer,globalcontainer_view);
      }
 
 
      template<typename M, typename GCView>
      void etadd (const M& localcontainer, GCView& globalcontainer_view) const
      {
 
        typedef typename M::value_type T;
 
        typedef typename GCView::RowIndexCache LFSVIndexCache;
        typedef typename GCView::ColIndexCache LFSUIndexCache;
 
        const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
        const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
 
        typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
        const LFSV& lfsv = lfsv_indices.localFunctionSpace();
 
        typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
        const LFSU& lfsu = lfsu_indices.localFunctionSpace();
 
        for (size_t i = 0; i<lfsv_indices.size(); ++i)
          for (size_t j = 0; j<lfsu_indices.size(); ++j)
            {
 
              if (localcontainer(lfsv,i,lfsu,j) == 0.0)
                continue;
 
              const bool constrained_v = lfsv_indices.isConstrained(i);
              const bool constrained_u = lfsu_indices.isConstrained(j);
 
              typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
              typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
 
              if (constrained_v)
                {
                  if (lfsv_indices.isDirichletConstraint(i))
                    continue;
 
                  for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
                    if (constrained_u)
                      {
                        if (lfsu_indices.isDirichletConstraint(j))
                          {
                            T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
                            if (value != 0.0)
                              globalcontainer_view.add(vcit->containerIndex(),j,value);
                          }
                        else
                          {
                            for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
                              {
                                T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
                                if (value != 0.0)
                                  globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),value);
                              }
                          }
                      }
                    else
                      {
                        T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
                        if (value != 0.0)
                          globalcontainer_view.add(vcit->containerIndex(),j,value);
                      }
                }
              else
                {
                  if (constrained_u)
                    {
                      if (lfsu_indices.isDirichletConstraint(j))
                        {
                          T value = localcontainer(lfsv,i,lfsu,j);
                          if (value != 0.0)
                            globalcontainer_view.add(i,j,value);
                        }
                      else
                        {
                          for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
                            {
                              T value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
                              if (value != 0.0)
                                globalcontainer_view.add(i,ucit->containerIndex(),value);
                            }
                        }
                    }
                  else
                    globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
                }
            }
      }
 
 
      template<typename Pattern, typename RI, typename CI>
      typename std::enable_if<
        std::is_same<RI,CI>::value
        >::type
      add_diagonal_entry(Pattern& pattern, const RI& ri, const CI& ci) const
      {
        if (ri == ci)
          pattern.add_link(ri,ci);
      }
 
      template<typename Pattern, typename RI, typename CI>
      typename std::enable_if<
        !std::is_same<RI,CI>::value
        >::type
      add_diagonal_entry(Pattern& pattern, const RI& ri, const CI& ci) const
      {}
 
      template<typename P, typename LFSVIndices, typename LFSUIndices, typename Index>
      void add_entry(P & globalpattern,
                     const LFSVIndices& lfsv_indices, Index i,
                     const LFSUIndices& lfsu_indices, Index j) const
      {
        typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
        typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
 
        const bool constrained_v = lfsv_indices.isConstrained(i);
        const bool constrained_u = lfsu_indices.isConstrained(j);
 
        add_diagonal_entry(globalpattern,
                           lfsv_indices.containerIndex(i),
                           lfsu_indices.containerIndex(j)
                           );
 
        if(!constrained_v)
          {
            if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
              {
                globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
              }
            else
              {
                for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
                  globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
              }
          }
        else
          {
            if (lfsv_indices.isDirichletConstraint(i))
              {
                globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
              }
            else
              {
                for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
                  {
                    if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
                      {
                        globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
                      }
                    else
                      {
                        for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
                          globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
                      }
                  }
              }
          }
      }
 
      template<typename GFSV, typename GC, typename C>
      void set_trivial_rows(const GFSV& gfsv, GC& globalcontainer, const C& c) const
      {
        typedef typename C::const_iterator global_row_iterator;
        for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
          globalcontainer.clear_row(cit->first,1);
      }
 
      template<typename GFSV, typename GC>
      void set_trivial_rows(const GFSV& gfsv, GC& globalcontainer, const EmptyTransformation& c) const
      {
      }
 
      template<typename GFSV, typename GC>
      void handle_dirichlet_constraints(const GFSV& gfsv, GC& globalcontainer) const
      {
        globalcontainer.flush();
        set_trivial_rows(gfsv,globalcontainer,*pconstraintsv);
        globalcontainer.finalize();
      }
 
      /* constraints */
      const CU* pconstraintsu;
      const CV* pconstraintsv;
      static CU emptyconstraintsu;
      static CV emptyconstraintsv;
    };
 
    template<typename B, typename CU, typename CV>
    CU LocalAssemblerBase<B,CU,CV>::emptyconstraintsu;
    template<typename B, typename CU, typename CV>
    CV LocalAssemblerBase<B,CU,CV>::emptyconstraintsv;
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH