dirichletconstraints.hh

#ifndef DUNE_ACFEM_DIRICHLETCONSTRAINTS_HH
#define DUNE_ACFEM_DIRICHLETCONSTRAINTS_HH
 
#include <dune/fem/function/common/scalarproducts.hh>
 
#include "../../models/boundaryindicator.hh"
#include "../../functions/boundarysupportedfunction.hh"
#include "../../functions/gridfunctionexpression.hh"
#include "bulkblockconstraints.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class JacobianOperator, class BoundaryValues, class Indicator = EntireBoundaryIndicatorType >
    class DirichletConstraints;
 
    template<class JacobianOperator, class BoundaryValues, class Indicator>
    class LocalDirichletConstraints;
 
    template<class JacobianOperator, class BoundaryValues, class Indicator>
    struct DirichletBlockConstraintsTraits
      : public
    BlockConstraintsDefaultStorageTraits<
      DirichletConstraints<JacobianOperator, BoundaryValues, Indicator>,
      LocalDirichletConstraints<JacobianOperator, BoundaryValues, Indicator> >
    {};
 
    template<class JacobianOperator, class BoundaryValues>
    struct DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType>
      : public BlockConstraintsDefaultTraits<DirichletConstraints<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType>,
                                             LocalDirichletConstraints<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType> >
    {};
 
    template<class JacobianOperator, class BoundaryValues>
    class DirichletConstraints<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType>
      : public DifferentiableBlockConstraintsOperatorInterface<JacobianOperator,
                                    DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType> >
    {
     public:
      typedef DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType> TraitsType;
      typedef typename TraitsType::LocalOperatorType LocalOperatorType;
      typedef JacobianOperator JacobianOperatorType;
      typedef typename JacobianOperatorType::DomainFunctionType DiscreteFunctionType;
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef BoundaryValues BoundaryValuesFunctionType;
      typedef EmptyBoundaryIndicatorType IndicatorType;
 
      DirichletConstraints(const DiscreteFunctionSpaceType& space,
                           const BoundaryValuesFunctionType& boundaryValues,
                           const IndicatorType& indicator)
      {}
 
      void operator()(const DiscreteFunctionType& arg, DiscreteFunctionType& result) const {}
      void jacobian(const DiscreteFunctionType& u, JacobianOperatorType& jOp) const {}
 
      LocalOperatorType localOperator(const EntityType& entity) const
      {
        return LocalOperatorType();
      }
 
      LocalOperatorType localOperator() const
      {
        return LocalOperatorType();
      }
 
      void rebuild() const {}
      void constrain(DiscreteFunctionType& w) const {}
      void zeroConstrain(DiscreteFunctionType& w) const {}
      size_t size() const { return 0; }
    };
 
 
    template<class JacobianOperator, class BoundaryValues>
    class LocalDirichletConstraints<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType>
    {
     protected:
      typedef DirichletConstraints<JacobianOperator, BoundaryValues, EmptyBoundaryIndicatorType> OperatorType;
      friend OperatorType;
 
      LocalDirichletConstraints(const OperatorType& op) {}
      LocalDirichletConstraints() {}
     public:
      typedef typename OperatorType::EntityType EntityType;
      typedef typename OperatorType::DiscreteFunctionType::LocalFunctionType LocalFunctionType;
      typedef JacobianOperator JacobianOperatorType;
 
      void init(const EntityType& entity) {}
      void operator()(const LocalFunctionType& u, LocalFunctionType& w) const {}
      template<class LocalMatrix>
      void jacobian(const LocalFunctionType& u /* unused */, LocalMatrix& jOpLocal) const {}
      void constrain(LocalFunctionType& u) const {}
      void zeroConstrain(LocalFunctionType& u) const {}
      size_t numConstrained() { return 0; }
      bool unconstrained() { return true; }
      bool totallyConstrained() const  { return false; }
    };
 
    template<class JacobianOperator, class BoundaryValues, class Indicator>
    class DirichletConstraints
      : public BulkBlockConstraints<JacobianOperator,
                                    DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, Indicator> >
    {
      typedef DirichletConstraints ThisType;
     public:
      typedef DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, Indicator> TraitsType;
     private:
      typedef BulkBlockConstraints<JacobianOperator, TraitsType> BaseType;
      typedef DifferentiableBlockConstraintsOperatorInterface<JacobianOperator, ThisType> InterfaceType;
     public:
      typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
      typedef typename BaseType::JacobianOperatorType JacobianOperatorType;
      typedef typename BaseType::LocalOperatorType LocalOperatorType;
      typedef typename BaseType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename BaseType::GridPartType GridPartType;
      typedef typename BaseType::EntityType EntityType;
 
      typedef BoundaryValues BoundaryValuesFunctionType;
      typedef decltype(*std::declval<Indicator>()) IndicatorType;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      // types for boundary treatment
      // ----------------------------
      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;
      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::localOperator;
 
      //using BaseType::operator();
      void operator()(const DiscreteFunctionType& arg, DiscreteFunctionType& result) const {
        BaseType::operator()(arg, result);
      }
 
      void rebuild() const
      {
        // std::cerr << "rebuild called "
        //           << sequence_
        //           << " "
        //           << space_.sequence()
        //           << std::endl;
        if (sequence_ != space_.sequence()) {
 
          BaseType::rebuildSlaveDofs();
 
          updateBoundaryValues(boundaryValues_(), affineOffset_);
 
          sequence_ = space_.sequence();
        }
      }
 
      void rebuildValues()
      {
        if (sequence_ != space_.sequence()) {
          BaseType::rebuildSlaveDofs();
        }
 
        updateBoundaryValues(boundaryValues_(), affineOffset_);
 
        sequence_ = space_.sequence();
      }
 
     protected:
      enum { localBlockSize = DiscreteFunctionSpaceType::localBlockSize };
 
      void updateBoundaryValues(const BoundaryValuesFunctionType& u, DiscreteFunctionType& w) const
      {
        const auto& mapper = space_.blockMapper();
        const int maxNumDofs = mapper.maxNumDofs();
 
        // resize flag vector with number of blocks and reset flags
        constrainedDofsMask_.resize(mapper.size());
        std::fill(constrainedDofsMask_.begin(), constrainedDofsMask_.end(), 0U);
 
        // only start search if Dirichlet boundary is present
        if (ExpressionTraits<IndicatorType>::isZero) {
          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 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;
              }
              ++constrainedDofsMask_[blockDofIndices[localBlock]];
              ++numConstrained_;
 
              int localDof = localBlock * localBlockSize;
              for(int l = 0; l < localBlockSize; ++l, ++localDof) {
                wLocal[localDof] = uValues[localDof];
              }
            }
          }
        }
      }
 
     protected:
      using BaseType::space_;
      using BaseType::affineOffset_;
      using BaseType::constrainedDofsMask_;
      using BaseType::sequence_;
      using BaseType::numConstrained_;
 
      const ExpressionStorage<BoundaryValuesFunctionType> boundaryValues_;
      const ExpressionStorage<IndicatorType> indicator_;
    };
 
    template<class JacobianOperator, class BoundaryValues, class Indicator>
    class LocalDirichletConstraints
      : public LocalBulkBlockConstraints<JacobianOperator,
                                         DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, Indicator> >
    {
      typedef DirichletBlockConstraintsTraits<JacobianOperator, BoundaryValues, Indicator> TraitsType;
      typedef LocalBulkBlockConstraints<JacobianOperator, TraitsType> BaseType;
     protected:
      typedef typename TraitsType::GlobalOperatorType OperatorType;
      typedef typename OperatorType::LocalObjectFactoryType LocalObjectFactoryType;
      friend LocalObjectFactoryType;
 
      LocalDirichletConstraints(const OperatorType& op)
        : BaseType(op)
      {}
     public:
      bool totallyConstrained() const
      {
        return false;
      }
    };
 
 
 
  } // ACFem::
 
} // Dune::
#endif