bulkblockconstraints.hh

#ifndef _DUNE_ACFEM_BULK_BLOCKCONSTRAINTS_HH_
#define _DUNE_ACFEM_BULK_BLOCKCONSTRAINTS_HH_
 
#include "blockconstraintsoperator.hh"
#include "../../common/functor.hh"
#include <dune/fem/misc/mpimanager.hh>
 
namespace Dune {
 
  namespace ACFem {
 
    template<class JacobianOperator, class Traits>
    class BulkBlockConstraints;
 
    template<class JacobianOperator, class Traits>
    class LocalBulkBlockConstraints;
 
    template<class JacobianOperator>
    struct BulkBlockConstraintsDefaultTraits
      : public
    BlockConstraintsDefaultStorageTraits<
      BulkBlockConstraints<JacobianOperator,
                           BulkBlockConstraintsDefaultTraits<JacobianOperator> >,
      LocalBulkBlockConstraints<JacobianOperator,
                                BulkBlockConstraintsDefaultTraits<JacobianOperator> > >
    {};
 
    template<class JacobianOperator, class Traits = BulkBlockConstraintsDefaultTraits<JacobianOperator> >
    class BulkBlockConstraints
      : public DifferentiableBlockConstraintsOperatorInterface<JacobianOperator, Traits>
    {
      typedef BulkBlockConstraints ThisType;
      typedef
      DifferentiableBlockConstraintsOperatorInterface<JacobianOperator, Traits>
      BaseType;
     public:
      typedef Traits TraitsType;
      typedef typename TraitsType::GlobalOperatorType GlobalOperatorType;
      typedef JacobianOperator JacobianOperatorType;
      typedef typename JacobianOperatorType::DomainFunctionType DiscreteFunctionType;
      typedef typename TraitsType::LocalOperatorType LocalOperatorType;
 
      typedef typename TraitsType::LocalObjectWrapperType LocalObjectWrapperType;
      typedef typename TraitsType::LocalObjectType LocalObjectType;
      typedef typename TraitsType::LocalObjectFactoryType LocalObjectFactoryType;
      typedef typename TraitsType::LocalObjectStorageType LocalObjectStorageType;
 
      typedef LocalBulkBlockConstraints<JacobianOperatorType, TraitsType> LocalBaseObjectType;
 
      friend LocalObjectType; // not necessarily identical to LocalBulkBlockConstraints
      friend LocalBaseObjectType;
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
     protected:
      typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
      typedef Dune::Fem::SlaveDofs<DiscreteFunctionSpaceType, BlockMapperType> SlaveDofsType;
      typedef typename SlaveDofsType::SingletonKey SlaveDofsKeyType;
      typedef Dune::Fem::SingletonList<SlaveDofsKeyType, SlaveDofsType> SlaveDofsProviderType;
      typedef typename BlockMapperType::GlobalKeyType GlobalKeyType;
 
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
     public:
 
      BulkBlockConstraints(const DiscreteFunctionType& values, const std::vector<unsigned>& mask)
        : space_(values.space()),
          localObjectFactory_(*static_cast<GlobalOperatorType *>(this)),
          localObjectStorage_(localObjectFactory_),
          affineOffset_(values),
          constrainedDofsMask_(mask),
          slaveDofMask_(space_.blockMapper().size(), false),
          numSlaveDofs_(space_.slaveDofs().size()-1),
          numConstrained_(std::count_if(constrainedDofsMask_.begin(), constrainedDofsMask_.end(),
                                        [](const unsigned& value) { return value > 0; })),
          sequence_(space_.sequence())
      {}
 
      BulkBlockConstraints(const DiscreteFunctionSpaceType& space)
        : space_(space),
          localObjectFactory_(*static_cast<GlobalOperatorType *>(this)),
          localObjectStorage_(localObjectFactory_),
          affineOffset_("Constraint Values", space_),
          constrainedDofsMask_(space_.blockMapper().size(), 0),
          slaveDofMask_(space_.blockMapper().size(), false),
          numSlaveDofs_(space_.slaveDofs().size()-1),
          numConstrained_(0),
          sequence_(-1)
      {}
 
      void rebuild() const
      {
        if (sequence_ != space_.sequence()) {
          rebuildSlaveDofs();
        }
 
        constrainedDofsMask_.resize(space_.blockMapper().size(), 0);
 
        numConstrained_ =
          std::count_if(constrainedDofsMask_.begin(), constrainedDofsMask_.end(),
                        [](const unsigned& value) { return value > 0; });
 
        sequence_ = space_.sequence();
      }
 
     protected:
      void rebuildSlaveDofs() const
      {
        const auto& slaveDofs = space_.slaveDofs();
        numSlaveDofs_ = slaveDofs.size() - 1;
 
        slaveDofMask_.resize(space_.blockMapper().size());
        slaveDofMask_.clear();
        for (unsigned i = 0; i < numSlaveDofs_; ++i) {
          slaveDofMask_[slaveDofs[i]] = true;
        }
      }
 
      bool isSlaveDof(GlobalKeyType index) const
      {
        assert(Fem::MPIManager::size() > 1 || numSlaveDofs_ == 0);
        return numSlaveDofs_ > 0 && slaveDofMask_[index];
      }
     public:
      void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &w) const
      {
        BaseType::rebuild(); // Update if out of sequence_ number
 
        if (size() == 0) {
          return;
        }
 
        auto uIt      = u.dbegin();
        auto wIt      = w.dbegin();
        auto offsetIt = affineOffset_.dbegin();
        auto maskIt   = constrainedDofsMask_.begin();
        auto maskEnd  = constrainedDofsMask_.end();
        while (maskIt != maskEnd) {
          if (*maskIt) {
            for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt, ++uIt) {
              *wIt = *uIt - *offsetIt;
            }
          } else {
            for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt, ++uIt);
          }
          ++maskIt;
        }
      }
 
      void jacobian(const DiscreteFunctionType& u, JacobianOperatorType& jOp) const
      {
        BaseType::rebuild(); // Update if out of sequence_ number
 
        if (size() == 0) {
          return;
        }
 
        const auto end = space_.end();
        for (auto it = space_.begin(); it != end; ++it) {
          auto localConstraints = localOperator(*it);
          auto jOpLocal = jOp.localMatrix(*it, *it);
 
        localConstraints.jacobian(u.localFunction(*it), jOpLocal);
        }
 
        jOp.communicate();
      }
 
      void constrain(DiscreteFunctionType& w) const
      {
        BaseType::rebuild(); // Update if out of sequence_ number
 
        if (size() == 0) {
          // nothing to do.
          return;
        }
 
        auto wIt      = w.dbegin();
        auto offsetIt = affineOffset_.dbegin();
        auto maskIt   = constrainedDofsMask_.begin();
        auto maskEnd  = constrainedDofsMask_.end();
        while (maskIt != maskEnd) {
          if (*maskIt) {
            for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt) {
              *wIt = *offsetIt;
            }
          } else {
            for (int l = 0; l < localBlockSize; ++l, ++wIt, ++offsetIt);
          }
          ++maskIt;
        }
      }
 
      void zeroConstrain(DiscreteFunctionType& w) const
      {
        BaseType::rebuild();
 
        if (size() == 0) {
          return;
        }
 
        auto wIt      = w.dbegin();
        auto maskIt   = constrainedDofsMask_.begin();
        auto maskEnd  = constrainedDofsMask_.end();
        while (maskIt != maskEnd) {
          if (*maskIt) {
            for (int l = 0; l < localBlockSize; ++l, ++wIt) {
              *wIt = 0;
            }
          } else {
            for (int l = 0; l < localBlockSize; ++l, ++wIt);
          }
          ++maskIt;
        }
      }
 
      size_t size() const {
        return numConstrained_;
      }
 
      LocalOperatorType localOperator(const EntityType& entity) const
      {
        return LocalOperatorType(entity, localObjectStorage_);
      }
 
      LocalOperatorType localOperator() const
      {
        return LocalOperatorType(localObjectStorage_);
      }
 
     protected:
      const DiscreteFunctionSpaceType& space_;
      mutable LocalObjectFactoryType localObjectFactory_;
      mutable LocalObjectStorageType localObjectStorage_;
      mutable DiscreteFunctionType affineOffset_;
      mutable std::vector<unsigned> constrainedDofsMask_;
      mutable std::vector<bool> slaveDofMask_;
      mutable std::size_t numSlaveDofs_;
      mutable unsigned numConstrained_;
      mutable int sequence_;
    };
 
    template<class JacobianOperator, class Traits>
    class LocalBulkBlockConstraints
    {
     protected:
      typedef Traits TraitsType;
      typedef BulkBlockConstraints<JacobianOperator, TraitsType> OperatorType;
      typedef JacobianOperator JacobianOperatorType;
      typedef typename JacobianOperatorType::RangeFunctionType DiscreteFunctionType;
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
      typedef typename BlockMapperType::GlobalKeyType GlobalKeyType;
 
      friend typename OperatorType::LocalObjectFactoryType;
 
      LocalBulkBlockConstraints(const OperatorType& op)
        : entity_(0),
          operator_(op),
          numberOfLocalBlocks_(0),
          constrainedDofs_(operator_.space_.blockMapper().maxNumDofs()),
          slaveDofs_(operator_.space_.blockMapper().maxNumDofs())
      {}
 
     public:
      typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
 
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
 
      void init(const EntityType& entity)
      {
        entity_ = &entity;
        numberOfLocalBlocks_ = operator_.space_.blockMapper().numDofs(*entity_);
 
#if 1
        // ... this way :::
        operator_.space_.blockMapper().mapEach(
          *entity_,
          makeTupleFunctor(ExtractFunctor<std::vector<unsigned>, std::vector<unsigned> >(constrainedDofs_, operator_.constrainedDofsMask_),
                           ExtractFunctor<std::vector<bool>, std::vector<bool> >(slaveDofs_, operator_.slaveDofMask_)));
#else
        // ... likewise ...
        operator_.space_.blockMapper().mapEach(
          *entity_,
          makePairFunctor(ExtractFunctor<std::vector<unsigned>, std::vector<unsigned> >(constrainedDofs_, operator_.constrainedDofsMask_),
                          ExtractFunctor<std::vector<bool>, std::vector<bool> >(slaveDofs_, operator_.slaveDofMask_)));
#endif
        numConstrained_ = std::count_if(constrainedDofs_.begin(), constrainedDofs_.end(),
                                        [](const unsigned& value) { return value > 0; });
      }
 
      bool unconstrained() const
      {
        return numConstrained_ == 0;
      }
 
      bool totallyConstrained() const
      {
        return numConstrained_ == numberOfLocalBlocks_;
      }
 
      size_t numConstrained() const
      {
        return numConstrained_;
      }
 
      void operator()(const LocalFunctionType& u, LocalFunctionType& w) const
      {
        if (unconstrained()) {
          return;
        }
 
        LocalFunctionType localOffset(operator_.affineOffset_.localFunction(*entity_));
 
        int localDof = 0;
        for (int localBlockDof = 0; localBlockDof < numberOfLocalBlocks_; ++localBlockDof) {
          if (constrainedDofs_[localBlockDof]) {
            if (slaveDofs_[localBlockDof]) {
              for (int l = 0; l < localBlockSize; ++l, ++localDof) {
                w[localDof] = 0;
              }
            } else {
              for (int l = 0; l < localBlockSize; ++l, ++localDof) {
                w[localDof] = u[localDof] - localOffset[localDof];
              }
            }
          } else {
            // Do nothing, operation is restricted to masked-out DOFs
            localDof += localBlockSize;
          }
        }
      }
 
      template<class LocalMatrix>
      void jacobian(const LocalFunctionType& u /* unused */, LocalMatrix& jOpLocal) const
      {
        if (unconstrained()) {
          return;
        }
 
        int localDof = 0;
        for (int localBlockDof = 0; localBlockDof < numberOfLocalBlocks_; ++localBlockDof) {
          unsigned constraintFactor = constrainedDofs_[localBlockDof]; // 0: unconstrained, > 0: number of elements sharing this DOF.
          if (constraintFactor > 0) {
            double value = slaveDofs_[localBlockDof] ? 0.0 : 1.0 / (double)constraintFactor;
            for (int l = 0; l < localBlockSize; ++l, ++localDof) {
              jOpLocal.clearRow(localDof); // indexed by non-block-DoF
              jOpLocal.set(localDof, localDof, value);
            }
          } else {
            localDof += localBlockSize;
          }
        }
      }
 
      void constrain(LocalFunctionType& w) const
      {
        if (unconstrained()) {
          return;
        }
 
        LocalFunctionType localOffset(operator_.affineOffset_.localFunction(*entity_));
 
        int localDof = 0;
        for (int localBlockDof = 0; localBlockDof < numberOfLocalBlocks_; ++localBlockDof) {
          if (constrainedDofs_[localBlockDof]) {
            if (slaveDofs_[localBlockDof]) {
              for (int l = 0; l < localBlockSize; ++l, ++localDof) {
                w[localDof] = 0;
              }
            } else {
              for (int l = 0; l < localBlockSize; ++l, ++localDof) {
                w[localDof] = localOffset[localDof];
              }
            }
          } else {
            // Do nothing, operation is restricted to masked-out DOFs
            localDof += localBlockSize;
          }
        }
      }
 
      void zeroConstrain(LocalFunctionType& w) const
      {
        if (unconstrained()) {
          return;
        }
 
        int localDof = 0;
        for (int localBlockDof = 0; localBlockDof < numberOfLocalBlocks_; ++localBlockDof) {
          if (constrainedDofs_[localBlockDof]) {
            for (int l = 0; l < localBlockSize; ++l, ++localDof) {
              w[localDof] = 0;
            }
          } else {
            // Do nothing, operation is restricted to masked-out DOFs
            localDof += localBlockSize;
          }
        }
      }
 
     protected:
      const EntityType* entity_;
      const OperatorType& operator_;
      int numberOfLocalBlocks_;
      std::vector<unsigned> constrainedDofs_;
      std::vector<bool> slaveDofs_;
      int numConstrained_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // _DUNE_ACFEM_EMPTY_BLOCKCONSTRAINTS_HH_