boundaryspace.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2013 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef BOUNDARY_SPACE_HH
#define BOUNDARY_SPACE_HH
 
#include <map>
#include <utility>
 
#include <dune/common/fvector.hh>
 
#include "fem/scalarspace.hh"
 
namespace Kaskade
{
  // forward declarations
  template <class,class> class ContinuousLagrangeMapper;
  template <class,class,class> class ContinuousLagrangeMapperImplementation;
  template <class,class> class DiscontinuousLagrangeMapper;
  template <class,class> class ContinuousHierarchicMapper;
  template <class,class> class DiscontinuousHierarchicMapper;
  template <class,class,class> class ContinuousHierarchicMapperImplementation;
  template <class,class> class ContinuousHierarchicExtensionMapper;
  template <class,class,class> class ContinuousHierarchicExtensionMapperImplementation;
 
  template <class,int,class> class LagrangeSimplexShapeFunctionSet;
  template <class,int,class> class HierarchicSimplexShapeFunctionSet;
  template <class,int,class> class HierarchicExtensionSimplexShapeFunctionSet;
 
  namespace BoundarySpace_Detail
  {
    using namespace ScalarSpaceDetail;
 
    template <class MapperImplementation> struct ChooseMapper;
 
    template <class Scalar, class GridView>
    struct ChooseMapper<ContinuousLagrangeMapper<Scalar,GridView> >
    {
      typedef ContinuousLagrangeMapperImplementation<Scalar,GridView,RestrictToBoundary> type;
    };
 
    struct BOUNDARY_MAPPERS_ARE_NOT_YET_IMPLEMENTED_FOR_DISCONTINUOUS_SPACES;
 
    template <class Scalar, class GridView>
    struct ChooseMapper<DiscontinuousLagrangeMapper<Scalar,GridView> >
    {
      typedef BOUNDARY_MAPPERS_ARE_NOT_YET_IMPLEMENTED_FOR_DISCONTINUOUS_SPACES type;
    };
 
    template <class Scalar, class GridView>
    struct ChooseMapper<ContinuousHierarchicMapper<Scalar,GridView> >
    {
      typedef ContinuousHierarchicMapperImplementation<Scalar,GridView,RestrictToBoundary> type;
    };
 
    template <class Scalar, class GridView>
    struct ChooseMapper<ContinuousHierarchicExtensionMapper<Scalar,GridView> >
    {
      typedef ContinuousHierarchicExtensionMapperImplementation<Scalar,GridView,RestrictToBoundary> type;
    };
 
    template <class Scalar, class GridView>
    struct ChooseMapper<DiscontinuousHierarchicMapper<Scalar,GridView> >
    {
      typedef BOUNDARY_MAPPERS_ARE_NOT_YET_IMPLEMENTED_FOR_DISCONTINUOUS_SPACES type;
    };
 
    template <class Mapper> struct GetSFDataType;
 
    template <class Scalar, class GridView>
    struct GetSFDataType<ContinuousLagrangeMapper<Scalar,GridView> >
    {
      typedef boost::compressed_pair<bool,Empty> type;
    };
 
    template <class Scalar, class GridView>
    struct GetSFDataType<ContinuousHierarchicMapper<Scalar,GridView> >
    {
      typedef boost::compressed_pair<bool,int> type;
    };
 
    template <class Scalar, class GridView>
    struct GetSFDataType<ContinuousHierarchicExtensionMapper<Scalar,GridView> >
    {
      typedef boost::compressed_pair<bool,int> type;
    };
 
    template <class Mapper> struct ShapeFunctionSetRestriction;
 
    template <class Scalar, class GridView, class Data>
    struct ShapeFunctionSetRestriction<ContinuousLagrangeMapperImplementation<Scalar,GridView,Data> >
    {
      typedef RestrictedShapeFunctionSet<LagrangeSimplexShapeFunctionSet<typename GridView::Grid::ctype,GridView::Grid::dimension,Scalar> > type;
    };
 
    template <class Scalar, class GridView, class Data>
    struct ShapeFunctionSetRestriction<ContinuousHierarchicMapperImplementation<Scalar,GridView,Data> >
    {
      typedef RestrictedShapeFunctionSet<HierarchicSimplexShapeFunctionSet<typename GridView::Grid::ctype,GridView::Grid::dimension,Scalar> > type;
    };
 
    template <class Scalar, class GridView, class Data>
    struct ShapeFunctionSetRestriction<ContinuousHierarchicExtensionMapperImplementation<Scalar,GridView,Data> >
    {
      typedef RestrictedShapeFunctionSet<HierarchicExtensionSimplexShapeFunctionSet<typename GridView::Grid::ctype,GridView::Grid::dimension,Scalar> > type;
    };
 
  }
 
  template <template <class,class> class Implementation_, class Scalar_, class GridView_>
  class [[deprecated]]
  BoundaryMapper : public UniformScalarMapper<typename BoundarySpace_Detail::ChooseMapper<Implementation_<Scalar_,GridView_>>::type,
                                              typename BoundarySpace_Detail::GetSFDataType<Implementation_<Scalar_,GridView_>>::type>
  {
    typedef typename BoundarySpace_Detail::ChooseMapper<Implementation_<Scalar_,GridView_>>::type Implementation;
    typedef typename BoundarySpace_Detail::GetSFDataType<Implementation_<Scalar_,GridView_>>::type DataType;
    typedef UniformScalarMapper<Implementation,DataType> Base;
 
  public:
    typedef Scalar_ Scalar;
    typedef GridView_ GridView;
    typedef typename BoundarySpace_Detail::ShapeFunctionSetRestriction<Implementation>::type ShapeFunctionSetImplementation;
    typedef int ConstructorArgument;
    static int const continuity = 0;
    static constexpr int dim = GridView::dimension;
 
    template <int m>
    struct Element
    {
      typedef FunctionSpaceElement<FEFunctionSpace<BoundaryMapper>,m> type;
    };
 
    BoundaryMapper(GridView const& gridView, int order):
      Base(Implementation(gridView,order))
    {
      assert(order >= 1);
    }
 
    BoundaryMapper(GridView const& gridView, int order, std::vector<int> const& boundaryIds_, std::vector<int> const& usedIds_)
     : Base(Implementation(gridView,order,ScalarSpaceDetail::RestrictToBoundary(boundaryIds_,usedIds_)))
    {
      assert(order>=1);
    }
 
    template <class T, class Functor>
    BoundaryMapper(GridView const& gridView, int order, std::vector<T> const& boundaryIds_, std::vector<T> const& usedIds_, Functor boundaryIdsToInts)
     : Base(Implementation(gridView,order,ScalarSpaceDetail::RestrictToBoundary(boundaryIdsToInts(boundaryIds_),boundaryIdsToInts(usedIds_))))
    {
      assert(order>=1);
    }
 
    bool inDomain(typename GridView::Intersection const& face) const
    {
      return ScalarSpaceDetail::considerFace(face,this->implementation.boundaryIds,this->implementation.usedIds);
    }
 
    int numFaces() {
      int nFaces = 0;
      forEachBoundaryFace(this->implementation.gridView(),[&](typename GridView::Intersection const& intersection) {
        if(inDomain(intersection)) ++nFaces;
      });
      return nFaces;
    }
  };
 
 
  namespace ScalarSpaceDetail
  {
    template <template <class,class,class> class MapperImplementation, class Scalar_, class GridView_, class SFData>
    class MapperPolicy<MapperImplementation<Scalar_,GridView_,RestrictToBoundary>, boost::compressed_pair<bool,SFData> >
    {
      typedef Scalar_ Scalar;
      typedef GridView_ GridView;
      typedef MapperImplementation<Scalar,GridView,RestrictToBoundary>                     Implementation;
      typedef typename Implementation::Grid                                                  Grid;
      typedef typename GridView::IndexSet                                                    IndexSet;
      typedef typename Grid::template Codim<0>::Entity                                      Cell;
      typedef typename Implementation::ShapeFunctionSet                                      ShapeFunctionSet;
      typedef boost::compressed_pair<bool,SFData>                                            Data;
      typedef typename BoundarySpace_Detail::ShapeFunctionSetRestriction<Implementation>::type  ShapeFunctionSet_Restricted;
 
    protected:
      static constexpr bool boundaryPolicy = true;
 
      explicit MapperPolicy(Implementation const& impl) : implementation(impl),
      insertionIndex(0), consideredCells(implementation.indexSet().size(0))
      {}
 
      template <class GlobalIndices, class SortedIndices>
      void initIndices(std::map<Dune::GeometryType,size_t>& startIndex, GlobalIndices& globIdx, SortedIndices& sortedIdx)
      {
        insertionIndex = 0;
        globIdx.clear();
        globIdx.resize(implementation.indexSet().size(0));
        sortedIdx.clear();
        sortedIdx.resize(implementation.indexSet().size(0));
 
        consideredCells.clear();
        consideredCells.resize(implementation.indexSet().size(0),false);
        shapeFunctionSets.clear();
        shapeFunctionSets.resize(implementation.indexSet().size(0));
 
        size_t n = 0;
        for (int codim=0; codim<=Grid::dimension; ++codim)
          for(auto const& geoType: implementation.indexSet().types(codim))
          {
            startIndex[geoType] = n;
            n += implementation.indexSet().size(geoType) * implementation.dofOnEntity(geoType);
          }
 
        indexMapper.clear();
        indexMapper.resize(n,std::make_pair(false,0));
 
        restrictCells();
      }
 
      template <class GlobalIndices, class SortedIndices>
      void computeIndices(Cell const& cell, std::map<Dune::GeometryType,size_t>& startIndex, GlobalIndices& globIdx, SortedIndices& sortedIdx, size_t cellIndex)
      {
        // No evaluation will be done on faces that do not intersect the boundary.
        // If evaluation should be allowed make sure to set a, possibly empty, vector using setRestriction()
        std::vector<int> idRestriction;
        ShapeFunctionSet& sf = static_cast<UniformScalarMapper<Implementation,Data>&>(*this).shapefunctions_non_const(cell);
        if(!consideredCells[cellIndex])
        {
          if(cell.type().isSimplex()) static_cast<ShapeFunctionSet_Restricted&>(sf).setRestriction(idRestriction);
          return;
        }
 
        int nominalOrder; // unused
 
        // step through all shape functions.
        Dune::GeometryType gt;
        int subentity, codim, indexInSubentity;
        size_t realLocalId = 0;
        for(size_t i=0; i<sf.size(); ++i)
        {
          std::tie(nominalOrder,codim,subentity,indexInSubentity) = sf[i].location();
 
          if(codim == 0) continue;
 
          Data data(false);
          assert(cell.type().isSimplex());
          implementation.entityIndex(cell,sf[i],i,gt,subentity,codim,indexInSubentity,data);
          if(data.first())
          {
            std::pair<bool,size_t> contains = std::make_pair(false,0);
            if(codim > 0)
            {
              // compute the global index of the subentity to which the shape function is associated
              int gt_idx = subIndex(implementation.indexSet(),cell,codim,subentity);
 
              std::map<Dune::GeometryType,size_t>::const_iterator mi = startIndex.find(gt);
              assert(mi!=startIndex.end());
 
              // On the boundary codim-0-entities do vanish
              // Shape functions associated with codim-1-entities are unique
              // thus we only have to check shape functions associated with entities
              // of codim > 1
              size_t unrestrictedIndex = mi->second+gt_idx*implementation.dofOnEntity(gt)+indexInSubentity;
              if(codim > 1) {
                contains = indexMapper[unrestrictedIndex];
                if(!contains.first) indexMapper[unrestrictedIndex] = std::make_pair(true,insertionIndex);
              } else {
                // for shapefunctions on codim 1 entities also store the restricted index (originally this map was only provided and used internally for codim >1 shapefunctions),
                indexMapper[unrestrictedIndex] = std::make_pair(true,insertionIndex);
              }
            }
 
            if(!contains.first) globIdx[cellIndex].push_back( boost::compressed_pair<size_t,Data>(insertionIndex,data) );
            else globIdx[cellIndex].push_back( boost::compressed_pair<size_t,Data>(contains.second,data) );
 
            sortedIdx[cellIndex].push_back( std::make_pair(globIdx[cellIndex].back().first(),realLocalId) );
 
            idRestriction.push_back(i);
 
            if(!contains.first) ++insertionIndex;
            ++realLocalId;
          }
        }
 
        if(cell.type().isSimplex()) static_cast<ShapeFunctionSet_Restricted&>(sf).setRestriction(idRestriction);
      }
 
      template <class Container, class ShapeFunctionSet>
      void initShapeFunctionSet(Container& sfs, ShapeFunctionSet const& sf, size_t cellIndex)
      {
        // TODO: probably it is not the most efficient idea to create for each cell an own restricted shape function set
        // the next line was replaced since it caused a memory leak
        //sfs[cellIndex] = new ShapeFunctionSet_Restricted(static_cast<ShapeFunctionSet_Restricted const&>(sf),nullptr); // cache shape function set
        assert(cellIndex < shapeFunctionSets.size());
        shapeFunctionSets[cellIndex] = std::make_unique<ShapeFunctionSet_Restricted>(static_cast<ShapeFunctionSet_Restricted const&>(sf),nullptr);
        sfs[cellIndex] = shapeFunctionSets[cellIndex].get();
      }
 
      Implementation implementation;
 
    public:
      std::pair<bool,size_t> unrestrictedToRestrictedIndex(std::vector<std::pair<bool,size_t>>::size_type unrestrictedIndex) const {
        assert(unrestrictedIndex < indexMapper.size());
        return indexMapper[unrestrictedIndex];
      }
 
    private:
      void restrictCells()
      {
        auto cend = implementation.gridView().template end<0>();
        for(auto cell = implementation.gridView().template begin<0>(); cell!=cend; ++cell)
        {
          size_t cellIndex = implementation.indexSet().index(*cell);
          auto fend = implementation.gridView().iend(*cell);
          for(auto face = implementation.gridView().ibegin(*cell); face!=fend; ++face)
            if(ScalarSpaceDetail::considerFace(*face,implementation.boundaryIds,implementation.usedIds))
              consideredCells[cellIndex] = true;
        }
      }
 
      size_t insertionIndex;
      std::vector<bool> consideredCells;
      // use the global numbering of standard spaces for fast identification of shared degrees of freedoms
      std::vector<std::pair<bool,size_t> > indexMapper;
      std::vector<std::unique_ptr<ShapeFunctionSet>> shapeFunctionSets; //store non-const copies of restricted shape function sets
    };
 
  } // end namespace ScalarSpaceDetail
} // end namespace Kaskade
 
#endif