axisalignedcubegeometry.hh

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:
 
#ifndef DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
#define DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
 
#include <bitset>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/diagonalmatrix.hh>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
 
namespace Dune {
 
  template <class CoordType, unsigned int dim, unsigned int coorddim>
  class AxisAlignedCubeGeometry
  {
 
 
  public:
 
    enum {mydimension = dim};
 
    enum {coorddimension = coorddim};
 
    typedef CoordType ctype;
 
    typedef FieldVector<ctype,dim> LocalCoordinate;
 
    typedef FieldVector<ctype,coorddim> GlobalCoordinate;
 
    typedef ctype Volume;
 
    typedef typename std::conditional<dim==coorddim,
        DiagonalMatrix<ctype,dim>,
        FieldMatrix<ctype,dim,coorddim> >::type JacobianTransposed;
 
    typedef typename std::conditional<dim==coorddim,
        DiagonalMatrix<ctype,dim>,
        FieldMatrix<ctype,coorddim,dim> >::type JacobianInverseTransposed;
 
    AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower,
                            const Dune::FieldVector<ctype,coorddim> upper)
      : lower_(lower),
        upper_(upper),
        axes_()
    {
      static_assert(dim==coorddim, "Use this constructor only if dim==coorddim!");
      // all 'true', but is never actually used
      axes_ = (1<<coorddim)-1;
    }
 
    AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower,
                            const Dune::FieldVector<ctype,coorddim> upper,
                            const std::bitset<coorddim>& axes)
      : lower_(lower),
        upper_(upper),
        axes_(axes)
    {
      assert(axes.count()==dim);
      for (size_t i=0; i<coorddim; i++)
        if (not axes_[i])
          upper_[i] = lower_[i];
    }
 
    AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower)
      : lower_(lower)
    {}
 
    GeometryType type() const
    {
      return GeometryTypes::cube(dim);
    }
 
    GlobalCoordinate global(const LocalCoordinate& local) const
    {
      GlobalCoordinate result;
      if (dim == coorddim) {        // fast case
        for (size_t i=0; i<coorddim; i++)
          result[i] = lower_[i] + local[i]*(upper_[i] - lower_[i]);
      } else if (dim == 0) {              // a vertex -- the other fast case
        result = lower_;          // hope for named-return-type-optimization
      } else {          // slow case
        size_t lc=0;
        for (size_t i=0; i<coorddim; i++)
          result[i] = (axes_[i])
                      ? lower_[i] + local[lc++]*(upper_[i] - lower_[i])
                      : lower_[i];
      }
      return result;
    }
 
    LocalCoordinate local(const GlobalCoordinate& global) const
    {
      LocalCoordinate result;
      if (dim == coorddim) {        // fast case
        for (size_t i=0; i<dim; i++)
          result[i] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
      } else if (dim != 0) {          // slow case
        size_t lc=0;
        for (size_t i=0; i<coorddim; i++)
          if (axes_[i])
            result[lc++] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
      }
      return result;
    }
 
    JacobianTransposed jacobianTransposed([[maybe_unused]] const LocalCoordinate& local) const
    {
      JacobianTransposed result;
 
      // Actually compute the result.  Uses different methods depending
      // on what kind of matrix JacobianTransposed is.
      jacobianTransposed(result);
 
      return result;
    }
 
    JacobianInverseTransposed jacobianInverseTransposed([[maybe_unused]] const LocalCoordinate& local) const
    {
      JacobianInverseTransposed result;
 
      // Actually compute the result.  Uses different methods depending
      // on what kind of matrix JacobianTransposed is.
      jacobianInverseTransposed(result);
 
      return result;
    }
 
    ctype integrationElement([[maybe_unused]] const LocalCoordinate& local) const
    {
      return volume();
    }
 
    GlobalCoordinate center() const
    {
      GlobalCoordinate result;
      if (dim==0)
        result = lower_;
      else {
        // Since lower_==upper_ for unused coordinates, this always does the right thing
        for (size_t i=0; i<coorddim; i++)
          result[i] = CoordType(0.5) * (lower_[i] + upper_[i]);
      }
      return result;
    }
 
    int corners() const
    {
      return 1<<dim;
    }
 
    GlobalCoordinate corner(int k) const
    {
      GlobalCoordinate result;
      if (dim == coorddim) {         // fast case
        for (size_t i=0; i<coorddim; i++)
          result[i] = (k & (1<<i)) ? upper_[i] : lower_[i];
      } else if (dim == 0) {         // vertex
        result = lower_;            // rely on named return-type optimization
      } else {                // slow case
        unsigned int mask = 1;
 
        for (size_t i=0; i<coorddim; i++) {
          if (not axes_[i])
            result[i] = lower_[i];
          else {
            result[i] = (k & mask) ? upper_[i] : lower_[i];
            mask = (mask<<1);
          }
        }
      }
 
 
      return result;
    }
 
    Volume volume() const
    {
      ctype vol = 1;
      if (dim == coorddim) {         // fast case
        for (size_t i=0; i<dim; i++)
          vol *= upper_[i] - lower_[i];
        // do nothing if dim == 0
      } else if (dim != 0) {         // slow case
        for (size_t i=0; i<coorddim; i++)
          if (axes_[i])
            vol *= upper_[i] - lower_[i];
      }
      return vol;
    }
 
    bool affine() const
    {
      return true;
    }
 
    friend Dune::Transitional::ReferenceElement< ctype, Dim<dim> > referenceElement ( const AxisAlignedCubeGeometry & /* geometry */ )
    {
      return ReferenceElements< ctype, dim >::cube();
    }
 
  private:
    // jacobianTransposed: fast case --> diagonal matrix
    void jacobianTransposed ( DiagonalMatrix<ctype,dim> &jacobianTransposed ) const
    {
      for (size_t i=0; i<dim; i++)
        jacobianTransposed.diagonal()[i] = upper_[i] - lower_[i];
    }
 
    // jacobianTransposed: slow case --> dense matrix
    void jacobianTransposed ( FieldMatrix<ctype,dim,coorddim> &jacobianTransposed ) const
    {
      if (dim==0)
          return;
 
      size_t lc = 0;
      for (size_t i=0; i<coorddim; i++)
        if (axes_[i])
          jacobianTransposed[lc++][i] = upper_[i] - lower_[i];
    }
 
    // jacobianInverseTransposed: fast case --> diagonal matrix
    void jacobianInverseTransposed ( DiagonalMatrix<ctype,dim> &jacobianInverseTransposed ) const
    {
      for (size_t i=0; i<dim; i++)
        jacobianInverseTransposed.diagonal()[i] = CoordType(1.0) / (upper_[i] - lower_[i]);
    }
 
    // jacobianInverseTransposed: slow case --> dense matrix
    void jacobianInverseTransposed ( FieldMatrix<ctype,coorddim,dim> &jacobianInverseTransposed ) const
    {
      if (dim==0)
          return;
 
      size_t lc = 0;
      for (size_t i=0; i<coorddim; i++)
        if (axes_[i])
          jacobianInverseTransposed[i][lc++] = CoordType(1.0) / (upper_[i] - lower_[i]);
    }
 
    Dune::FieldVector<ctype,coorddim> lower_;
 
    Dune::FieldVector<ctype,coorddim> upper_;
 
    std::bitset<coorddim> axes_;
  };
 
} // namespace Dune
#endif