simplex.cc

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
#define DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
 
// This file is part of DUNE, a Distributed and Unified Numerics Environment
// This file is copyright (C) 2005 Jorrit Fahlke <jorrit@jorrit.de>
// This file is licensed under version 2 of the GNU General Public License,
// with a special "runtime exception."  See COPYING at the top of the source
// tree for the full licence.
 
#include <algorithm>
 
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
 
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include "base.cc"
 
namespace Dune {
 
  namespace RefinementImp {
 
    namespace Simplex {
 
      // //////////////////
      //
      //
 
 
      [[deprecated("Use factorial from dune-common's math.hh")]]
      inline int factorial(int n)
      {
        int prod = 1;
        for(int i = 1; i <= n; ++i)
          prod *= i;
        return prod;
      }
 
      [[deprecated("Use binomial from dune-common's math.hh")]]
      inline int binomial(int upper, int lower)
      {
        lower = std::min( lower, upper - lower );
        if(lower < 0)
          return 0;
        int prod = 1;
        for(int i = upper - lower; i < upper; ++i)
          prod *= (i+1);
        return prod / Dune::factorial(lower);
      }
 
      template<int dimension>
      int pointIndex(const FieldVector<int, dimension> &point)
      {
        int index = 0;
        for(int i = 0; i < dimension; ++i)
          index += Dune::binomial(dimension-i + point[i]-1, dimension-i);
        return index;
      }
 
      template<int n>
      FieldVector<int, n> getPermutation(int m)
      {
        FieldVector<int, n> perm;
        for(int i = 0; i < n; ++i)
          perm[i] = i;
 
        int base = 1;
        for(int i = 1; i <= n; ++i)
          base *= i;
 
        for(int i = n; i > 0; --i) {
          base /= i;
          int d = m / base;
          m %= base;
          int t = perm[i-1]; perm[i-1] = perm[i-1-d]; perm[i-1-d] = t;
        }
        return perm;
      }
 
#if 0
      Has to be checked
      // calculate the index of a permutation
      template<int n>
      int getPermIndex(const FieldVector<int, n>& test) // O(n^2)
      {
        int m = 0;
        FieldVector<int, n> perm;
        for(int i = 0; i < n; ++i)
          perm[i] = i;
 
        int base = 1;
        for(int i = 1; i <= n; ++i)
          base *= i;
 
        for(int i = n; i > 0; --i) {
          base /= i;
          int d;
          for(d = 0; d < i; ++d)
            if(test[i-1] == perm[i-1-d])
              break;
          m += d * base;
          int d = m / base;
          m %= base;
          perm[i-1-d] = perm[i-1];
        }
      }
#endif
 
      // map between the reference simplex and some arbitrary kuhn simplex (denoted by it's permutation)
      template<int dimension, class CoordType>
      FieldVector<CoordType, dimension>
      referenceToKuhn( 
        FieldVector<CoordType, dimension> point,
        const FieldVector<int, dimension> &kuhn)
      {
        for(int i = dimension - 1; i > 0; --i)
          point[kuhn[i-1]] += point[kuhn[i]];
        return point;
      }
 
      template<int dimension, class CoordType>
      FieldVector<CoordType, dimension>
      kuhnToReference( 
        FieldVector<CoordType, dimension> point,
        const FieldVector<int, dimension> &kuhn)
      {
        for(int i = 0; i < dimension - 1; ++i)
          point[kuhn[i]] -= point[kuhn[i+1]];
        return point;
      }
 
 
 
      // /////////////////////////////////////////
      //
      // refinement implementation for simplices
      //
 
      template<int dimension_, class CoordType>
      class RefinementImp
      {
      public:
        constexpr static int dimension = dimension_;
        typedef CoordType ctype;
 
        template<int codimension>
        struct Codim;
        typedef typename Codim<dimension>::SubEntityIterator VertexIterator;
        typedef FieldVector<CoordType, dimension> CoordVector;
        typedef typename Codim<0>::SubEntityIterator ElementIterator;
        typedef FieldVector<int, dimension+1> IndexVector;
 
        static int nVertices(int nIntervals);
        static VertexIterator vBegin(int nIntervals);
        static VertexIterator vEnd(int nIntervals);
 
        static int nElements(int nIntervals);
        static ElementIterator eBegin(int nIntervals);
        static ElementIterator eEnd(int nIntervals);
      };
 
      template<int dimension, class CoordType>
      template<int codimension>
      struct RefinementImp<dimension, CoordType>::Codim
      {
        class SubEntityIterator;
        // We don't need the caching, but the uncached MultiLinearGeometry has bug FS#1209
        typedef Dune::CachedMultiLinearGeometry<CoordType,dimension-codimension,dimension> Geometry;
      };
 
      template<int dimension, class CoordType>
      int
      RefinementImp<dimension, CoordType>::
      nVertices(int nIntervals)
      {
        return Dune::binomial(dimension + nIntervals, (int)dimension);
      }
 
      template<int dimension, class CoordType>
      typename RefinementImp<dimension, CoordType>::VertexIterator
      RefinementImp<dimension, CoordType>::
      vBegin(int nIntervals)
      {
        return VertexIterator(nIntervals);
      }
 
      template<int dimension, class CoordType>
      typename RefinementImp<dimension, CoordType>::VertexIterator
      RefinementImp<dimension, CoordType>::
      vEnd(int nIntervals)
      {
        return VertexIterator(nIntervals, true);
      }
 
      template<int dimension, class CoordType>
      int
      RefinementImp<dimension, CoordType>::
      nElements(int nIntervals)
      {
        return Dune::power(nIntervals, int(dimension));
      }
 
      template<int dimension, class CoordType>
      typename RefinementImp<dimension, CoordType>::ElementIterator
      RefinementImp<dimension, CoordType>::
      eBegin(int nIntervals)
      {
        return ElementIterator(nIntervals);
      }
 
      template<int dimension, class CoordType>
      typename RefinementImp<dimension, CoordType>::ElementIterator
      RefinementImp<dimension, CoordType>::
      eEnd(int nIntervals)
      {
        return ElementIterator(nIntervals, true);
      }
 
      // //////////////
      //
      // The iterator
      //
 
      template<int dimension, class CoordType, int codimension>
      class RefinementIteratorSpecial;
 
      // vertices
 
      template<int dimension, class CoordType>
      class RefinementIteratorSpecial<dimension, CoordType, dimension>
      {
      public:
        typedef RefinementImp<dimension, CoordType> Refinement;
        typedef typename Refinement::CoordVector CoordVector;
        typedef typename Refinement::template Codim<dimension>::Geometry Geometry;
        typedef RefinementIteratorSpecial<dimension, CoordType, dimension> This;
 
        RefinementIteratorSpecial(int nIntervals, bool end = false);
 
        void increment();
        bool equals(const This &other) const;
 
        CoordVector coords() const;
        Geometry geometry () const;
 
        int index() const;
      protected:
        typedef FieldVector<int, dimension> Vertex;
 
        int size;
        Vertex vertex;
      };
 
      template<int dimension, class CoordType>
      RefinementIteratorSpecial<dimension, CoordType, dimension>::
      RefinementIteratorSpecial(int nIntervals, bool end)
        : size(nIntervals)
      {
        vertex[0] = (end) ? size + 1 : 0;
        for(int i = 1; i < dimension; ++ i)
          vertex[i] = 0;
      }
 
      template<int dimension, class CoordType>
      void
      RefinementIteratorSpecial<dimension, CoordType, dimension>::
      increment()
      {
        assert(vertex[0] <= size);
        for(int i = dimension - 1; i >= 0; --i) {
          ++vertex[i];
          if(i == 0 || vertex[i] <= vertex[i-1])
            break;
          else
            vertex[i] = 0;
        }
      }
 
      template<int dimension, class CoordType>
      bool
      RefinementIteratorSpecial<dimension, CoordType, dimension>::
      equals(const This &other) const
      {
        return size == other.size && vertex == other.vertex;
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, dimension>::CoordVector
      RefinementIteratorSpecial<dimension, CoordType, dimension>::
      coords() const
      {
        Vertex ref = kuhnToReference(vertex, getPermutation<dimension>(0));
 
        CoordVector coords;
        for(int i = 0; i < dimension; ++i)
          coords[i] = CoordType(ref[i]) / size;
        return coords;
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, dimension>::Geometry
      RefinementIteratorSpecial<dimension, CoordType, dimension>::geometry () const
      {
        std::vector<CoordVector> corners(1);
        corners[0] = (CoordVector)vertex;
        return Geometry(GeometryTypes::vertex, corners);
      }
 
      template<int dimension, class CoordType>
      int
      RefinementIteratorSpecial<dimension, CoordType, dimension>::
      index() const
      {
        return pointIndex(vertex);
      }
 
      // elements
 
      template<int dimension, class CoordType>
      class RefinementIteratorSpecial<dimension, CoordType, 0>
      {
      public:
        typedef RefinementImp<dimension, CoordType> Refinement;
        typedef typename Refinement::IndexVector IndexVector;
        typedef typename Refinement::CoordVector CoordVector;
        typedef typename Refinement::template Codim<0>::Geometry Geometry;
        typedef RefinementIteratorSpecial<dimension, CoordType, 0> This;
 
        RefinementIteratorSpecial(int nIntervals, bool end = false);
 
        void increment();
        bool equals(const This &other) const;
 
        IndexVector vertexIndices() const;
        int index() const;
        CoordVector coords() const;
 
        Geometry geometry () const;
 
      private:
        CoordVector global(const CoordVector &local) const;
 
      protected:
        typedef FieldVector<int, dimension> Vertex;
        constexpr static int nKuhnIntervals = Dune::factorial(dimension);
 
        Vertex origin;
        int kuhnIndex;
        int size;
        int index_;
      };
 
      template<int dimension, class CoordType>
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      RefinementIteratorSpecial(int nIntervals, bool end)
        : kuhnIndex(0), size(nIntervals), index_(0)
      {
        for(int i = 0; i < dimension; ++i)
          origin[i] = 0;
        if(end) {
          index_ = Refinement::nElements(nIntervals);
          origin[0] = size;
        }
      }
 
      template<int dimension, class CoordType>
      void
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      increment()
      {
        assert(origin[0] < size);
 
        ++index_;
 
        while(1) {
          ++kuhnIndex;
          if(kuhnIndex == nKuhnIntervals) {
            kuhnIndex = 0;
            // increment origin
            for(int i = dimension - 1; i >= 0; --i) {
              ++origin[i];
              if(i == 0 || origin[i] <= origin[i-1])
                break;
              else
                origin[i] = 0;
            }
          }
 
          // test whether the current simplex has any corner outside the kuhn0 simplex
          FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
          Vertex corner = origin;
          bool outside = false;
          for(int i = 0; i < dimension; ++i) {
            // next corner
            ++corner[perm[i]];
            if(perm[i] > 0)
              if(corner[perm[i]] > corner[perm[i]-1]) {
                outside = true;
                break;
              }
          }
          if(!outside)
            return;
        }
      }
 
      template<int dimension, class CoordType>
      bool
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      equals(const This &other) const
      {
        return size == other.size && index_ == other.index_;
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, 0>::IndexVector
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      vertexIndices() const
      {
        IndexVector indices;
        FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
        Vertex vertex = origin;
        indices[0] = pointIndex(vertex);
        for(int i = 0; i < dimension; ++i) {
          ++vertex[perm[i]];
          indices[i+1] = pointIndex(vertex);
        }
        if (kuhnIndex%2 == 1)
          for(int i = 0; i < (dimension+1)/2; ++i) {
            int t = indices[i];
            indices[i] = indices[dimension-i];
            indices[dimension-i] = t;
          }
        return indices;
      }
 
      template<int dimension, class CoordType>
      int
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      index() const
      {
        return index_;
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, 0>::CoordVector
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      coords() const
      {
        return global(ReferenceElements<CoordType, dimension>
                      ::simplex().position(0,0));
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, 0>::Geometry
      RefinementIteratorSpecial<dimension, CoordType, 0>::geometry () const
      {
        std::vector<CoordVector> corners(dimension+1);
        auto refelem = ReferenceElements<CoordType, dimension>::simplex();
        for(int i = 0; i <= dimension; ++i)
          corners[i] = global(refelem.position(i, dimension));
        return Geometry(refelem.type(), corners);
      }
 
      template<int dimension, class CoordType>
      typename RefinementIteratorSpecial<dimension, CoordType, 0>::CoordVector
      RefinementIteratorSpecial<dimension, CoordType, 0>::
      global(const CoordVector &local) const {
        CoordVector v =
          referenceToKuhn(local, getPermutation<dimension>(kuhnIndex));
        v += origin;
        v /= (typename CoordVector::value_type)size;
        return kuhnToReference(v, getPermutation<dimension>(0));
      }
 
      // common
 
      template<int dimension, class CoordType>
      template<int codimension>
      class RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator
        : public ForwardIteratorFacade<typename RefinementImp<dimension, CoordType>::template Codim<codimension>::SubEntityIterator, int>,
          public RefinementIteratorSpecial<dimension, CoordType, codimension>
      {
      public:
        typedef RefinementImp<dimension, CoordType> Refinement;
 
        SubEntityIterator(int nIntervals, bool end = false);
      };
 
#ifndef DOXYGEN
 
      template<int dimension, class CoordType>
      template<int codimension>
      RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator::
      SubEntityIterator(int nIntervals, bool end)
        : RefinementIteratorSpecial<dimension, CoordType, codimension>(nIntervals, end)
      {}
 
#endif
 
    } // namespace Simplex
 
  } // namespace RefinementImp
 
 
  namespace RefinementImp {
 
    // ///////////////////////
    //
    // The refinement traits
    //
 
#ifndef DOXYGEN
    template<unsigned topologyId, class CoordType, unsigned coerceToId,
        int dim>
    struct Traits<
        topologyId, CoordType, coerceToId, dim,
        typename std::enable_if<
            ((GeometryTypes::simplex(dim).id() >> 1) ==
             (topologyId >> 1) &&
             (GeometryTypes::simplex(dim).id() >> 1) ==
             (coerceToId >> 1)
            )>::type
        >
    {
      typedef Simplex::RefinementImp<dim, CoordType> Imp;
    };
#endif
 
 
  } // namespace RefinementImp
 
} // namespace Dune
 
#endif //DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC