dirichletconstraints.hh

#ifndef DUNE_ACFEM_DIRICHLETCONSTRAINTS_HH
#define DUNE_ACFEM_DIRICHLETCONSTRAINTS_HH
 
#include <dune/fem/function/common/scalarproducts.hh>
 
#include "../../models/indicators/boundaryindicator.hh"
#include "../../models/modeltraits.hh"
#include "../../functions/functions.hh"
#include "bulkblockconstraints.hh"
#include "../../models/expressions.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace, class BoundaryValues, class Indicator = void, class SFINAE = void>
    class DirichletConstraints;
 
    namespace {
 
      template<class BoundaryValues, class SFINAE = void>
      struct HasEmptySupport
        : FalseType
      {};
 
      template<class BoundaryValues>
      struct HasEmptySupport<
        BoundaryValues,
        std::enable_if_t<(IsPDEModel<BoundaryValues>::value
                          && !EffectiveModelTraits<BoundaryValues>::template Exists<PDEModel::dirichlet>::value)> >
        : TrueType
      {};
 
    }
 
    template<class DiscreteFunctionSpace, class Model>
    class DirichletConstraints<DiscreteFunctionSpace, Model,
                               void,
                               std::enable_if_t<true
                                                && IsPDEModel<Model>::value
                                                && !HasEmptySupport<Model>::value> >
      : public BulkBlockConstraints<DiscreteFunctionSpace>
    {
      typedef DirichletConstraints ThisType;
      typedef BulkBlockConstraints<DiscreteFunctionSpace> BaseType;
     public:
      typedef typename BaseType::ConstraintsStorageType ConstraintsStorageType;
      typedef typename BaseType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename BaseType::GridPartType GridPartType;
      typedef typename BaseType::EntityType EntityType;
 
      using ModelType = EffectiveModel<Model>;
      using ModelTraits = EffectiveModelTraits<ModelType>;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      // types for boundary treatment
      // ----------------------------
      typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
 
      //typedef typename BlockMapperType::GlobalKeyType GlobalKeyType;
      typedef std::size_t GlobalKeyType;
 
      DirichletConstraints(const DiscreteFunctionSpaceType& space,
                           const ModelType& model)
        : BaseType(space), model_(model)
      {}
 
      using BaseType::size;
      using BaseType::jacobian;
      using BaseType::constrain;
      using BaseType::zeroConstrain;
      using BaseType::communicate;
 
      //using BaseType::operator();
      //      template<class DomainDiscreteFunction, class RangeDiscreteFunction>
      //void operator()(const DomainDiscreteFunction& arg, RangeDiscreteFunctionType& result) const {
      //BaseType::operator()(arg, result);
      //}
 
      void update() const
      {
        if (sequence_ == space_.sequence()) {
          return;
        }
 
        updateBoundaryValues();
 
        communicate();
 
        sequence_ = space_.sequence();
      }
 
     protected:
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
 
      void updateBoundaryValues() const
      {
        const auto& mapper = space_.blockMapper();
        const int maxNumDofs = mapper.maxNumDofs();
 
        // resize flag vector with number of blocks and reset flags
        mask_.resize(mapper.size());
        std::fill(mask_.begin(), mask_.end(), 0U);
        values_.clear();
 
        // only start search if Dirichlet boundary is present
        if (!ModelTraits::template Exists<PDEModel::dirichlet>::value) {
          numConstrained_ = 0;
          return;
        }
 
        // pre-allocate local objects
        std::vector<GlobalKeyType> blockDofIndices;
        std::vector<bool> blockDofFaceMask;
        blockDofIndices.reserve(maxNumDofs);
        blockDofFaceMask.reserve(maxNumDofs);
        auto uValues = temporaryLocalFunction(space_);
        auto model   = model_;
        auto uLocal  = modelDirichletValues(model, space_.gridPart());
        auto wLocal  = values_.localFunction();
 
        numConstrained_ = 0;
        for (const auto& entity : space_) {
          for (const auto& intersection : intersections(space_.gridPart(), entity)) {
 
            // if intersection is with boundary, adjust data
            if (intersection.neighbor() || !intersection.boundary()) {
              continue;
            }
 
            uLocal.bind(entity); // need to bind to  intersection.inside()
            const auto indicator = uLocal.classifyIntersection(intersection);
            const auto& localMask = indicator.second;
 
            if (!indicator.first || localMask.none()) {
              continue;
            }
 
            const auto numConstrained = localMask.count();
 
            // actual number of DoFs
            const int numDofs = mapper.numDofs(entity);
 
            // init local stuff
            blockDofIndices.resize(numDofs, GlobalKeyType(-1));
            blockDofFaceMask.resize(numDofs, false);
            uValues.init(entity);
            wLocal.init(entity);
 
            // get local -> global mapping
            space_.blockMapper().map(entity, blockDofIndices);
 
            // get mask of affected (block)-DoFs
            // get face number of boundary intersection
            const int face = intersection.indexInInside();
            mapper.onSubEntity(entity, face, 1 /* codim */, blockDofFaceMask);
 
            // call "natural" interpolation for _ALL_ DoFs
            space_.interpolation(entity)(uLocal, uValues);
 
            for (int localBlock = 0; localBlock < numDofs; ++localBlock) {
              if (!blockDofFaceMask[localBlock]) {
                continue;
              }
              auto& blockMask = mask_[blockDofIndices[localBlock]];
              blockMask |= localMask;
              auto numBlockConstrained = blockMask.count();
              numConstrained_ += 2*numConstrained - numBlockConstrained;
              int localDof = localBlock * localBlockSize;
              for(int l = 0; l < localBlockSize; ++l, ++localDof) {
                wLocal[localDof] = uValues[localDof];
              }
            } // block-dof loop
          } // intersection loop
        } // mesh-loop
      }  // method
 
     protected:
      using BaseType::space_;
      using BaseType::values_;
      using BaseType::mask_;
      using BaseType::sequence_;
      using BaseType::numConstrained_;
      ModelType model_;
    };
 
    template<class DiscreteFunctionSpace, class BoundaryValues, class Indicator>
    class DirichletConstraints<DiscreteFunctionSpace, BoundaryValues, Indicator,
                               std::enable_if_t<(IsWrappableByConstLocalFunction<BoundaryValues>::value
                                                 && IsBoundaryIndicator<Indicator>::value
                                                 && !IndicatorTraits<Indicator>::emptySupport
                                 )> >
      : public BulkBlockConstraints<DiscreteFunctionSpace>
    {
      typedef DirichletConstraints ThisType;
      typedef BulkBlockConstraints<DiscreteFunctionSpace> BaseType;
     public:
      typedef typename BaseType::ConstraintsStorageType ConstraintsStorageType;
      typedef typename BaseType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename BaseType::GridPartType GridPartType;
      typedef typename BaseType::EntityType EntityType;
 
      using IndicatorType = Indicator;
      typedef BoundaryValues BoundaryValuesFunctionType;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      // types for boundary treatment
      // ----------------------------
      typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
 
      //typedef typename BlockMapperType::GlobalKeyType GlobalKeyType;
      typedef std::size_t GlobalKeyType;
 
      DirichletConstraints(const DiscreteFunctionSpaceType& space,
                           const BoundaryValuesFunctionType& boundaryValues,
                           const IndicatorType& indicator)
        : BaseType(space)
        , boundaryValues_(boundaryValues)
        , indicator_(indicator)
      {}
 
      using BaseType::size;
      using BaseType::jacobian;
      using BaseType::constrain;
      using BaseType::zeroConstrain;
      using BaseType::communicate;
 
      //using BaseType::operator();
      //      template<class DomainDiscreteFunction, class RangeDiscreteFunction>
      //void operator()(const DomainDiscreteFunction& arg, RangeDiscreteFunctionType& result) const {
      //BaseType::operator()(arg, result);
      //}
 
      void update() const
      {
        if (sequence_ == space_.sequence()) {
          return;
        }
 
        updateBoundaryValues(boundaryValues_(), values_);
 
        communicate();
 
        sequence_ = space_.sequence();
      }
 
     protected:
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
 
      void updateBoundaryValues(const BoundaryValuesFunctionType& u, ConstraintsStorageType& w) const
      {
        const auto& mapper = space_.blockMapper();
        const int maxNumDofs = mapper.maxNumDofs();
 
        // resize flag vector with number of blocks and reset flags
        mask_.resize(mapper.size());
        std::fill(mask_.begin(), mask_.end(), 0U);
 
        // pre-allocate local objects
        std::vector<GlobalKeyType> blockDofIndices;
        std::vector<bool> blockDofFaceMask;
        blockDofIndices.reserve(maxNumDofs);
        blockDofFaceMask.reserve(maxNumDofs);
        auto uValues = temporaryLocalFunction(space_);
        auto uLocal = localFunction(u);
        auto wLocal = localFunction(w);
 
        numConstrained_ = 0;
        for (const auto& entity : space_) {
          for (const auto& intersection : intersections(space_.gridPart(), entity)) {
 
            // if intersection is with boundary, adjust data
            if (intersection.neighbor() || !intersection.boundary()) {
              continue;
            }
 
            // Check whether this segment belongs to a Dirichlet boundary ...
            if (!indicator_.applies(intersection)) {
              continue; // Nope
            }
 
            // actual number of DoFs
            const int numDofs = mapper.numDofs(entity);
 
            // init local stuff
            blockDofIndices.resize(numDofs, GlobalKeyType(-1));
            blockDofFaceMask.resize(numDofs, false);
            uValues.init(entity);
            uLocal.init(entity, intersection);
            wLocal.init(entity);
 
            // get local -> global mapping
            space_.blockMapper().map(entity, blockDofIndices);
 
            // get mask of affected (block)-DoFs
            // get face number of boundary intersection
            const int face = intersection.indexInInside();
            mapper.onSubEntity(entity, face, 1 /* codim */, blockDofFaceMask);
 
            // call "natural" interpolation for _ALL_ DoFs
            space_.interpolation(entity)(uLocal, uValues);
 
            for (int localBlock = 0; localBlock < numDofs; ++localBlock) {
              if (!blockDofFaceMask[localBlock]) {
                continue;
              }
              if (!mask_[blockDofIndices[localBlock]].any()) {
                // FIXME: comonent wise
                mask_[blockDofIndices[localBlock]].set();
                numConstrained_ += mask_[blockDofIndices[localBlock]].count();
              }
              int localDof = localBlock * localBlockSize;
              for(int l = 0; l < localBlockSize; ++l, ++localDof) {
                wLocal[localDof] = uValues[localDof];
              }
            }
          } // intersection loop
        } // mesh-loop
      }  // method
 
     protected:
      using BaseType::space_;
      using BaseType::values_;
      using BaseType::mask_;
      using BaseType::sequence_;
      using BaseType::numConstrained_;
 
      const BoundaryValuesFunctionType boundaryValues_;
      const IndicatorType indicator_;
    };
 
    template<class DiscreteFunctionSpace, class BoundaryValues, class Indicator>
    class DirichletConstraints<DiscreteFunctionSpace, BoundaryValues, Indicator,
                               std::enable_if_t<(HasEmptySupport<BoundaryValues>::value
                                                 || IndicatorTraits<Indicator>::emptySupport
      )> >
    {
     public:
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef BoundaryValues BoundaryValuesFunctionType;
      typedef EmptyBoundaryIndicator IndicatorType;
 
      DirichletConstraints(const DiscreteFunctionSpaceType& space,
                           const BoundaryValuesFunctionType& boundaryValues)
      {}
 
      template<class DomainDiscreteFunction, class RangeDiscreteFunction>
      void operator()(const DomainDiscreteFunction& arg, RangeDiscreteFunction& result) const {}
 
      template<class DomainDiscreteFunction, class JacobianOperator>
      void jacobian(const DomainDiscreteFunction& u, JacobianOperator& jOp) const {}
 
      template<class DomainDiscreteFunction, class JacobianOperator>
      void applyToOperator(const DomainDiscreteFunction& u, JacobianOperator& jOp) const
      {}
 
      template<class DomainDiscreteFunction>
      void constrain(DomainDiscreteFunction& w) const {}
 
      template<class DomainDiscreteFunction>
      void zeroConstrain(DomainDiscreteFunction& w) const {}
 
      size_t size() const { return 0; }
      void reset() const {}
      void update() const {}
      void updateValues() const {}
    };
 
 
 
  } // ACFem::
 
} // Dune::
#endif