bulkblockconstraints.hh

#ifndef _DUNE_ACFEM_BULK_BLOCKCONSTRAINTS_HH_
#define _DUNE_ACFEM_BULK_BLOCKCONSTRAINTS_HH_
 
#include "../../common/functor.hh"
#include "../../common/dofmappertupledatahandle.hh"
#include <dune/fem/operator/common/differentiableoperator.hh>
#include <dune/fem/function/blockvectorfunction/blockvectorfunction.hh>
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace,
             template<class T, class Alloc = std::allocator<T> > class MaskStorage = std::vector>
    class BulkBlockConstraints
    {
      typedef BulkBlockConstraints ThisType;
     public:
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      // AdaptiveDiscreteFunction does not work with our tuple data-handle.
      //typedef Fem::AdaptiveDiscreteFunction<DiscreteFunctionSpaceType> ConstraintsStorageType;
      typedef Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpaceType> ConstraintsStorageType;
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef std::bitset<localBlockSize> MaskType;
      typedef MaskStorage<MaskType> MaskStorageType;
     protected:
      typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
      typedef typename BlockMapperType::GlobalKeyType GlobalKeyType;
 
     public:
 
      BulkBlockConstraints(const DiscreteFunctionSpaceType& space)
        : space_(space),
          values_("Constraint Values", space_),
          mask_(space_.blockMapper().size(), 0),
          numConstrained_(0),
          sequence_(-1)
      {}
 
      virtual ~BulkBlockConstraints()
      {}
 
      virtual void update() const = 0;
#if 0
      {
        if (!force && sequence_ == space_.sequence()) {
          return;
        }
 
        mask_.resize(space_.blockMapper().size(), 0);
 
        numConstrained_ =
          std::count_if(mask_.begin(), mask_.end(),
                        [](const MaskType& value) { return value.count(); });
 
        communicate();
 
        sequence_ = space_.sequence();
      }
#endif
 
      virtual void updateValues() const {
        reset();
        update();
      }
 
      template<class DomainDiscreteFunction, class RangeDiscreteFunction>
      void operator()(const DomainDiscreteFunction& u, RangeDiscreteFunction& w) const
      {
        checkUpdate();
 
        if (size() == 0) {
          return;
        }
 
        auto uIt      = u.dbegin();
        auto wIt      = w.dbegin();
        auto offsetIt = values_.dbegin();
        for (const auto& mask : mask_) {
          for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt, ++uIt) {
            if (mask[l]) {
              *wIt = *uIt - *offsetIt;
            }
          }
        }
      }
 
      template<class DomainDiscreteFunction, class JacobianOperator>
      void jacobian(const DomainDiscreteFunction& u, JacobianOperator& jOp) const
      {
        checkUpdate();
 
        if (size() == 0) {
          return;
        }
 
        std::vector<size_t> rows;
        rows.reserve(size());
        size_t r = 0;
        for (const auto& mask : mask_) {
          for (int l = 0; l < localBlockSize; ++l) {
            if (mask[l]) {
              rows.push_back(r);
            }
            ++r;
          }
        }
        jOp.setUnitRows(rows);
      }
 
      template<class DomainDiscreteFunction, class JacobianOperator>
      void applyToOperator(const DomainDiscreteFunction& u, JacobianOperator& jOp) const
      {
        jacobian(u, jOp);
      }
 
      template<class DomainDiscreteFunction>
      void constrain(DomainDiscreteFunction& w) const
      {
        checkUpdate();
 
        if (size() == 0) {
          // nothing to do.
          return;
        }
 
        auto wIt      = w.dbegin();
        auto offsetIt = values_.dbegin();
        for (const auto& mask : mask_) {
          for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt) {
            if (mask[l]) {
              *wIt = *offsetIt;
            }
          }
        }
      }
 
      template<class DomainDiscreteFunction>
      void zeroConstrain(DomainDiscreteFunction& w) const
      {
        checkUpdate();
 
        if (size() == 0) {
          return;
        }
 
        auto wIt = w.dbegin();
        for (const auto& mask : mask_) {
          for (int l = 0; l < localBlockSize; ++l, ++wIt) {
            if (mask[l]) {
              *wIt = 0;
            }
          }
        }
      }
 
      size_t size() const {
        return numConstrained_;
      }
 
      void reset() const {
        sequence_ = -1;
      }
 
     protected:
      void checkUpdate() const
      {
        if (space_.sequence() != sequence_) {
          update();
        }
      }
 
      virtual void communicate() const
      {
        if (space_.gridPart().comm().size() <= 1) {
          return; // serial run
        }
 
        typedef
          DofMapperTupleDataHandleTraits<MaskStorageType,
                                         typename ConstraintsStorageType::DofVectorType>
          CommTraitsType;
        typedef typename CommTraitsType::DataItemType DataItemType;
        typedef typename CommTraitsType::DataItemScatterType DataItemScatterType;
 
        auto dataHandle = dofMapperTupleDataHandle(
          space_,
          [](const DataItemType& a, const DataItemScatterType& b) {
            for (int l = 0; l < localBlockSize; ++l) {
              // copy over data from a if it is constraint. Otherwise
              // keep our data. If both bits are set check for
              // consistency, otherwise we end up with
              // "constraint-ping-pong" with inconsistent values.
              assert(std::get<0>(b)[l] && std::get<0>(a)[l] &&
                     std::abs(std::get<1>(b)[l] - std::get<1>(a)[l]) < 1e-8);
              if (std::get<0>(a)[l]) {
                std::get<1>(b)[l] = std::get<1>(a)[l];
              }
            }
            std::get<0>(b) |= std::get<0>(a);
          },
          mask_,
          values_.dofVector());
 
        space_.gridPart().communicate(dataHandle,
                                      DiscreteFunctionSpaceType::GridPartType::indexSetInterfaceType,
                                      Dune::ForwardCommunication);
      }
 
     protected:
      const DiscreteFunctionSpaceType& space_;
      mutable ConstraintsStorageType values_;
      mutable MaskStorageType mask_;
      mutable unsigned numConstrained_;
      mutable int sequence_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // _DUNE_ACFEM_EMPTY_BLOCKCONSTRAINTS_HH_