partitionedspace.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) 2023-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef PARTITIONEDSPACE_HH
#define PARTITIONEDSPACE_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/lagrangespace.hh"
#include "fem/views.hh"
 
namespace Kaskade
{
  // forward declarations
  template <class,class> class ContinuousLagrangeBoundaryMapper;
 
  namespace PartitionedSpaceDetail
  {
    // an empty class for assotiating no additional data with localized ansatz functions
    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();
      }
    };
 
    // --------------------------------------------------------------------------------------------
 
    // Dummy tagger assigning the same tag to all cells, leading to a globally continuous space.
    struct SingleTagger
    {
      template <class Cell>
      int operator()(Cell const& cell) const
      {
        auto x = cell.geometry().center();
        return x[0]+x[1] < 1? 0: 1;
      }
    };
    
    // --------------------------------------------------------------------------------------------
 
   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 PartitionedSpaceDetail
  
  // --------------------------------------------------------------------------------------------
  // --------------------------------------------------------------------------------------------
 
  template <class SFData>
  using UniformPartitionedMapperCombinerArgument = RangeView<typename std::vector<boost::compressed_pair<size_t,SFData>>::const_iterator>;
 
  template <class Implementation,
            class Tagger = PartitionedSpaceDetail::SingleTagger,
            class SFData = PartitionedSpaceDetail::Empty>
  class UniformPartitionedMapper
  {
  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:
    // For each shape function (or localized ansatz function) on a cell we store both
    // the global index and (optionally) some additional data. For hierarchic mappers, e.g.,
    // this might be some sign for edge functions. Since often this is not needed, and
    // SFData is empty, we use the compressed pair memory optimization to avoid overhead.
    typedef boost::compressed_pair<size_t,SFData> Data;
 
    // To extract the global index from the Data compressed pair.
    typedef PartitionedSpaceDetail::CompressedFirst<Data> First;
 
    // An iterator extracting the global indices of localized ansatz functions.
    typedef boost::transform_iterator<First,typename std::vector<Data>::const_iterator> GlobalIndexIterator;
 
    // An iterator for accessing (global,local) index pairs of localized ansatz functions.
    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;
 
    UniformPartitionedMapper(Implementation const& impl, Tagger const& tagger_)
    : tagger(tagger_)
    , 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 nDofs;
    }
 
    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);
    }
 
    // TODO: Check where I'm used. doc me or remove me.
    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];    
    }
 
    // TODO: Check where I'm used. doc me or remove me.
    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();
 
      // First allocate the  memory needed for storing global indices
      // to prevent frequent reallocation.
      size_t const nCells = indexSet.size(0);
      sfs.resize(nCells);
      globIdx.resize(nCells);
      sortedIdx.resize(nCells);
 
 
      // In count, for any tag and any geometry type that occurs in the indexSet
      // we store a (sorted) list of the global entity indices included.
      // We request that there is a fixed number of dofs associated with each
      // tag-subentity type incidence, stored as second entry,
      // such that equal length index blocks can be assigned
      // to each such incidence. The later layout of dofs is thus as
      // follows (example continuous 2D simplex mesh with three triangles (two with tag 0),
      // order 3):
      // tag 0: [Interior Triangle Nodes, Interior Edge Nodes, Vertex Nodes]
      // tag 1: [Interior Triangle Nodes, Interior Edge Nodes, Vertex Nodes]
      // with substructuring
      // [ [[t1],[t2]], [[e1a,e1b],[e2a,e2b],[e3a,e3b],[e4a,e4b],[e5a,e5b]], [v1,v2,v3,v4],
      // [[t3]], [[e5a,e5b],[e6a,e6b],[e7a,e7b]], [v3,v4,v5] ]
      // where the first row corresponds to tag 0 and the second row to tag 1.
      // The start index for the dofs on the tag-type incidence is stored along
      // as third entry.
      using Tag = decltype(tagger(Cell()));
      using Key = std::pair<Tag,Dune::GeometryType>;
      using IncidenceData = std::tuple<std::vector<size_t>,int,size_t>;
      std::map<Key,IncidenceData> count;
 
 
      // In a two-stage approach, we first count the number of dofs for each tag-type incidence
      // by stepping through all cells, then compute the start indices for each such incidence.
      order = 0;
      for(auto const& cell : elements(gridView))
      {
        size_t const cellIndex = indexSet.index(cell);
        auto tag = tagger(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.
        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);
 
          // Find the the data for (tag,geometry type) incidences, inserting an empty datum
          // if there is none so far.
          auto mi = count.insert(std::pair(Key(tag,gt),
                                           IncidenceData{{},implementation.dofOnEntity(gt),0})).first;
          assert(mi!=count.end());
 
          // register the concrete entity (possibly several times)
          std::get<0>(mi->second).push_back(gt_idx);
        }
      }
 
      // Compute start indices by building partial sums. For that we need to
      // have the number of entities in each tag-type incidence, which we obtain
      // by making the list of actual entity indices unique.
      size_t start = 0;
      for (auto& entry: count)
      {
        // TODO: sorting might be expensive. Do profiling.
        auto& idxs = std::get<0>(entry.second);
        int const dofPerEntity = std::get<1>(entry.second);
 
        std::sort(idxs.begin(),idxs.end());
        idxs.erase(std::unique(idxs.begin(),idxs.end()),idxs.end());
 
        std::get<2>(entry.second) = start;   // enter the start position
        start += dofPerEntity*idxs.size();   // add the count to obtain next start position
      }
 
      // Now that we have seen all, start contains the total number of dofs.
      nDofs = start;
 
      // 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);
        auto tag = tagger(cell);
 
        // Skip this cell if there are no dofs associated.
        if (globIdx[cellIndex].empty())
          continue;
 
        // For each (local) shape function, have it's global index computed and noted.
        auto const& localSfs = *sfs[cellIndex];
        for(size_t i=0; i<localSfs.size(); ++i)
        {
          Dune::GeometryType gt;
          int subentity, codim, indexInSubentity;
          SFData sfData;
          implementation.entityIndex(cell,localSfs[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);
 
          // Find the count of (tag,geometry type) incidences, inserting zero
          // if there is none so far.
          auto mi = count.find(Key(tag,gt));
          assert(mi!=count.end());
 
          auto const& [eIdx,dofPerEntity,start] = mi->second;
          auto pos = std::lower_bound(begin(eIdx),end(eIdx),gt_idx);
          assert(pos != end(eIdx));
          int globalIndex = start + (pos-begin(eIdx))*dofPerEntity;
 
          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                              nDofs = 0;
 
    // 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;
 
    // The provider of cell tags.
    Tagger                              tagger;
 
  protected:
    Implementation                      implementation;
  };
 
 
  // --------------------------------------------------------------------------------------------
  // --------------------------------------------------------------------------------------------
 
  template <class Tagger, class GV, class ScalarType=double>
  class PiecewiseContinuousLagrangeMapper
  : public UniformPartitionedMapper<ContinuousLagrangeMapperImplementation<ScalarType,GV>,Tagger>
  {
    typedef ContinuousLagrangeMapperImplementation<ScalarType,GV> Implementation;
    typedef UniformPartitionedMapper<Implementation,Tagger> Base;
 
  public:
    typedef ScalarType Scalar;
    typedef GV GridView;
    typedef typename Base::ShapeFunctionSet ShapeFunctionSet;
    typedef int ConstructorArgument;
    static int const continuity = -1;
 
    template <int m>
    using Element_t = FunctionSpaceElement<FEFunctionSpace<PiecewiseContinuousLagrangeMapper>,m>;
 
    template <int m>
    struct [[deprecated]] Element
    {
      typedef FunctionSpaceElement<FEFunctionSpace<PiecewiseContinuousLagrangeMapper>,m> type;
    };
 
    PiecewiseContinuousLagrangeMapper(GridView const& gridView, int order, Tagger const& tagger)
    : Base(Implementation(gridView,order),tagger)
    {
      assert(order > 0);
    }
  };
 
} // end of namespace Kaskade
 
#endif