scalarspace.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-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SCALARSPACE_HH
#define SCALARSPACE_HH
 
#include <algorithm>
#include <functional>
#include <map>
#include <vector>
 
#include <boost/compressed_pair.hpp>
#include <boost/iterator/transform_iterator.hpp>
 
#include "fem/firstless.hh"
#include "fem/gridBasics.hh"
#include "fem/gridcombinatorics.hh"
#include "fem/views.hh"
 
namespace Kaskade
{
  // forward declarations
  template <class,class> class ContinuousLagrangeBoundaryMapper;
 
  namespace ScalarSpaceDetail
  {
    struct Empty {};
 
    template <class Pair>
    struct CompressedFirst
    {
      typedef typename Pair::first_const_reference result_type;
 
      typename Pair::first_const_reference operator()(Pair const& pair) const
      {
        return pair.first();
      }
    };
 
    // Overload here as compressed pair's entries are accessed by method call. 
    struct FirstLess
    {
      template <class Data>
      bool operator()(boost::compressed_pair<size_t,Data> const& a, boost::compressed_pair<size_t,Data> const& b)
      {
        return a.first() < b.first();
      }
    };
    
   inline bool onSimplexFace(int dim, int codim, int id, int localFaceId)
    {
      assert(dim<=3);
 
      if(codim==0) return false;                      // cell is never part of a face
      if(codim==1) return id==localFaceId;            // ... face only if it is the same face
 
      if(dim==2)
      {
        if(localFaceId==0) return (id==0 || id==1);
        if(localFaceId==1) return (id==0 || id==2);
        if(localFaceId==2) return (id==2 || id==1);
      }
      else if(dim==3)
      {
        // edges
        if(codim==2)
        {
          if(localFaceId==0) return (id>=0 && id<3);
          if(localFaceId==1) return (id==0 || id==3 || id==4);
          if(localFaceId==2) return (id==1 || id==3 || id==5);
          if(localFaceId==3) return (id==2 || id==4 || id==5);
        }
        // vertices
        if(codim==3)
        {
          if(localFaceId==0) return (id>=0 && id<3);
          if(localFaceId==1) return (id==0 || id==1 || id==3);
          if(localFaceId==2) return (id==0 || id==2 || id==3);
          if(localFaceId==3) return (id==1 || id==2 || id==3);
        }
        return false;
      }
 
      return false;                                   // never get here
    }
 
    static constexpr int defaultIndex = -94279;
 
    template <class Face>
    int getBoundaryId(Face const& face, std::vector<int> const& boundaryIds)
    {
      size_t boundarySegmentIndex = face.boundarySegmentIndex();
      if(boundarySegmentIndex < boundaryIds.size()) return boundaryIds[boundarySegmentIndex];
      return defaultIndex;
    }
 
    inline bool usedId(int id, std::vector<int> const& usedIds)
    {
      if (id==defaultIndex) return false;
      return std::find(usedIds.begin(), usedIds.end(), id) != usedIds.end();
    }
 
    template <class Face>
    inline bool considerFace(Face const& face, std::vector<int> const& boundaryIds, std::vector<int> const& usedIds)
    {
      if(boundaryIds.empty()) return face.boundary();
      return (face.boundary() && usedId(getBoundaryId(face,boundaryIds),usedIds));
    }
 
    struct RestrictToBoundary
    {
      RestrictToBoundary() : boundaryIds(), usedIds() {}
      RestrictToBoundary(RestrictToBoundary const& other) : boundaryIds(other.boundaryIds), usedIds(other.usedIds) {}
      RestrictToBoundary(std::vector<int> const& boundaryIds_, std::vector<int> const& usedIds_) : boundaryIds(boundaryIds_), usedIds(usedIds_)
      {}
 
      template <class Data, class GridView, class Cell>
      void treatBoundary(Data& data, GridView const& gridView, Cell const& cell, int codim, int subentity) const
      {
        // ignore shape functions that are associated with codim-0 entities
        assert(cell.type().isSimplex());
        if(codim > 0)
          for(auto const& face: intersections(gridView,cell))
            // only consider boundary faces
            if(considerFace(face,boundaryIds,usedIds)
               && onSimplexFace(GridView::dimension,codim,subentity,face.indexInInside()))
                data.first() = true;
      }
 
      std::vector<int> const boundaryIds;
      std::vector<int> const usedIds;
    };
 
    template <class Policy>
    constexpr bool isRestriction()
    {
      return std::is_same<Policy,RestrictToBoundary>::value;
    }
 
    struct AllShapeFunctions
    {
      template <class Data, class GridView, class Cell> void treatBoundary(Data&, GridView const&, Cell const&, int, int) const {}
    };
  } // End of namespace ScalarSpaceDetail
  
  // --------------------------------------------------------------------------------------------
  // --------------------------------------------------------------------------------------------
 
  template <class SFData>
  using UniformScalarMapperCombinerArgument = RangeView<typename std::vector<boost::compressed_pair<size_t,SFData>>::const_iterator>;
 
  template <class Implementation, class SFData = ScalarSpaceDetail::Empty>
  class UniformScalarMapper
  {
  public:
    typedef typename Implementation::Grid              Grid;
    using Cell = Kaskade::Cell<Grid>;
    typedef typename Implementation::ShapeFunctionSet  ShapeFunctionSet;
    typedef typename Implementation::Converter         Converter;
    typedef typename Implementation::Combiner          Combiner;
    typedef typename Implementation::Scalar            Scalar;
    typedef typename Implementation::GridView          GridView;
    typedef typename GridView::IndexSet                IndexSet;
    typedef std::pair<size_t,int>                      IndexPair;
 
  private:
    typedef boost::compressed_pair<size_t,SFData> Data;
 
    typedef ScalarSpaceDetail::CompressedFirst<Data> First;
 
    typedef boost::transform_iterator<First,typename std::vector<Data>::const_iterator> GlobalIndexIterator;
    typedef std::vector<IndexPair>::const_iterator                                      SortedIndexIterator;
 
    static constexpr int dim = Grid::dimension;
 
  public:
    typedef RangeView<GlobalIndexIterator> GlobalIndexRange;
    typedef RangeView<SortedIndexIterator> SortedIndexRange;
 
    static bool const globalSupport = false;
 
    UniformScalarMapper(Implementation const& impl) : implementation(impl)
    {
      update();
    }
 
    GridView const& gridView() const
    {
      return implementation.gridView();
    }
 
    int maxOrder() const
    {
      return order;
    }
    
    
    GlobalIndexRange initGlobalIndexRange() const
    {
      return GlobalIndexRange(GlobalIndexIterator(globIdx[0].begin(),First()),
                              GlobalIndexIterator(globIdx[0].begin(),First()));
    }
 
 
    GlobalIndexRange globalIndices(Cell const& cell) const
    {
      return globalIndices(implementation.indexSet().index(cell));
    }
 
    GlobalIndexRange globalIndices(size_t n) const
    {
      return GlobalIndexRange(GlobalIndexIterator(globIdx[n].begin(),First()),
                              GlobalIndexIterator(globIdx[n].end(),  First()));
    }
 
    static SortedIndexRange initSortedIndexRange()
    {
      static std::vector<IndexPair> dummy; // empty
      return SortedIndexRange(dummy.begin(),dummy.end());
    }
 
    SortedIndexRange sortedIndices(Cell const& cell) const
    {
      return sortedIndices(implementation.indexSet().index(cell));
    }
 
    SortedIndexRange sortedIndices(size_t n) const
    {
      return SortedIndexRange(sortedIdx[n].begin(),sortedIdx[n].end());
    }
 
    size_t size() const { return n; }
 
    ShapeFunctionSet const& shapefunctions(Cell const& cell, bool contained=false) const
    {
      if (contained || implementation.indexSet().contains(cell))
        return shapefunctions(implementation.indexSet().index(cell));
      else
        return implementation.shapeFunctions(cell);
    }
 
    ShapeFunctionSet& shapefunctions_non_const(Cell const& cell)
    {
      return shapefunctions_non_const(implementation.indexSet().index(cell));
    }
 
    ShapeFunctionSet const& shapefunctions(size_t n) const
    {
      return *sfs[n];    
    }
 
    ShapeFunctionSet& shapefunctions_non_const(size_t n)
    {
      return *sfs[n];
    }
 
    ShapeFunctionSet const& lowerShapeFunctions(Cell const& cell) const
    {
      if (globalIndices(cell).empty())
        return implementation.emptyShapeFunctionSet();
      else
        return implementation.lowerShapeFunctions(cell);
    }
 
    Combiner combiner(Cell const& cell, size_t index) const
    {
      assert(implementation.indexSet().index(cell)==index);
      return Combiner(rangeView(globIdx[index].begin(),globIdx[index].end()),
                      shapefunctions(index));
    }
 
 
    // just returns the index, since nothing is restricted
    std::pair<bool,size_t> unrestrictedToRestrictedIndex(size_t unrestrictedIndex)
    {
      return std::make_pair(true, unrestrictedIndex);
    }
 
    void update()
    {
      GridView const& gridView = implementation.gridView();
      IndexSet const& indexSet = implementation.indexSet();
 
      // Notify the implementation of the update request.
      implementation.update();
 
      // For each codimension (i.e. type of subentity) compute the
      // number of global ansatz functions as well as an accumulated
      // index into an array of all ansatz functions.
      startIndex.clear();
 
      size_t ndofs = 0;
      for (int codim=0; codim<=Grid::dimension; ++codim)
        for(auto const& geoType: implementation.indexSet().types(codim))
        {
          startIndex[geoType] = ndofs;
 
          size_t numOfEntity = implementation.size(geoType);
          size_t  dofOnEntity = implementation.dofOnEntity(geoType);
          ndofs += numOfEntity * dofOnEntity;
        }
 
      // Store the total number of degrees of freedom.
      n = ndofs;
 
      sfs.resize(indexSet.size(0));
      order = 0;
 
      // Precompute and cache all the global indices. First allocate the
      // memory needed to prevent frequent reallocation.
      globIdx.resize(implementation.indexSet().size(0));
      sortedIdx.resize(implementation.indexSet().size(0));
      
      // Step through all cells and compute the global ansatz function
      // indices of all shape functions on that cell.
      for(auto const& cell : elements(gridView))
      {
        size_t const cellIndex = indexSet.index(cell);
 
        // Store a shape function set even if the current cell does not belong to our
        // support subdomain. This is necessary because we return a *reference*
        // to a shape function set for all cells.
        ShapeFunctionSet const& sf = implementation.shapeFunctions(cell);
        sfs[cellIndex] = &sf;
        
 
        // Skip this cell and let it's set of dofs empty if it is not in our support.
        if (!implementation.inSupport(cell))
        {
          globIdx[cellIndex].clear();     // We must clear the indices explicitly, as on refinement
          sortedIdx[cellIndex].clear();   // there may be leftovers from the previous grid.
          continue;
        }
 
        // Track the maximum shape function order encountered on any cell.
        order = std::max(order,sf.order());
 
        size_t localNumberOfShapeFunctions = sf.size();
        globIdx[cellIndex].resize(localNumberOfShapeFunctions);
        sortedIdx[cellIndex].resize(localNumberOfShapeFunctions);
 
        // For each (local) shape function, have it's global index computed and noted.
        // TODO: pull common parts out of this loop
        for(size_t i=0; i<localNumberOfShapeFunctions; ++i)
        {
          Dune::GeometryType gt;
          int subentity, codim, indexInSubentity;
          SFData sfData;
          implementation.entityIndex(cell,sf[i],i,gt,subentity,codim,indexInSubentity,sfData);
 
          // compute the global index of the subentity to which the shape function is associated
          int gt_idx = implementation.subIndex(cell,subentity,codim);
 
          auto mi = startIndex.find(gt);
          assert(mi!=startIndex.end());
 
          int globalIndex = mi->second+gt_idx*implementation.dofOnEntity(gt) + indexInSubentity;
 
          globIdx[cellIndex][i] =  Data(globalIndex,sfData);
          sortedIdx[cellIndex][i] = std::make_pair(globIdx[cellIndex][i].first(),i);
        }
 
 
        // sort the index pairs according to the global index
        std::sort(sortedIdx[cellIndex].begin(),sortedIdx[cellIndex].end(),FirstLess());
 
        // make sure that any assigned shape function forms at most one ansatz function.
        // This is a functionality restriction, but seems quite reasonable as a debugging check.
        // Note that not all shape functions need not be assigned. If a shape function is mapped     //Jakob: Are negative indices for unused shape function really allowed? Are'nt they simply not present in globIdx, such that the check idx[j]<0 below is nonsense?
        // to negative global index, it takes not part at all (e.g. in hp-methods).
#ifndef NDEBUG
        std::vector<size_t> idx(globIdx[cellIndex].size());
        for (int j=0; j<idx.size(); ++j)
          idx[j] = globIdx[cellIndex][j].first();
        std::sort(idx.begin(),idx.end());
        for (int j=1; j<idx.size(); ++j)
          assert(idx[j]>idx[j-1] || idx[j]<0);
#endif
      }
    }
 
 
  private:
    // The number of dofs.
    size_t                              n;
 
    // In startIndex, for any geometry type that occurs in the indexSet
    // there is an index stored, such that the dofs associated with
    // nodes on subentities of this type start at this index. We request
    // that there is a fixed number of dofs associated with each
    // subentity type, such that equal length index blocks are assigned
    // to each subentity of a given type. The layout of dofs is thus as
    // follows (example continuous 2D simplex mesh with two triangles,
    // order 4):
    //
    // [Interior Triangle Nodes, Interior Edge Nodes, Vertex Nodes] with substructuring
    // [ [[t1a,t1b,t1c],[t2a,t2b,t2c]], [[e1a,e1b],[e2a,e2b],[e3a,e3b],[e4a,e4b],[e5a,e5b]], [v1,v2,v3,v4] ]
    //     | start index triangles        | start index edges                                 | start index vertices
    // localDofs(triangles)=3, localDofs(edges)=2, localDofs(vertices)=1
    std::map<Dune::GeometryType,size_t> startIndex;
 
    // In globIdx, for each cell there are the global ansatz function indices on that cell.
    std::vector<std::vector<Data>>      globIdx;
 
    // In sortedIdx, for each cell there is a vector of (global index, local index) pairs, sorted ascendingly
    // by the global index. This could be computed outside on demand, but the sorting takes time and may be
    // amortized over multiple assembly passes/matrix blocks if cached here.
    std::vector<std::vector<IndexPair>> sortedIdx;
 
    // In sfs, for each cell there is a pointer to the shape function set on this cell cached.
    std::vector<ShapeFunctionSet const*> sfs;
 
    // The maximum polynomial ansatz order of shape functions on the cells.
    int                                 order;
 
 
  protected:
    Implementation                      implementation;
  };
} // end of namespace Kaskade
 
#endif