mmgeometry.hh 5.28 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1 2 3 4 5 6 7 8 9 10 11
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MULTIMESH_LOCALGEOMETRY_HH
#define DUNE_MULTIMESH_LOCALGEOMETRY_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
#include <dune/grid/common/geometry.hh>

namespace Dune
{
12
  template <int mydim, int coorddim, class HostGrid>
Praetorius, Simon's avatar
Praetorius, Simon committed
13
  class MultiMeshLocalGeometry
14
      : public GeometryDefaultImplementation<mydim, coorddim, HostGrid, MultiMeshLocalGeometry>
Praetorius, Simon's avatar
Praetorius, Simon committed
15 16
  {
  private:
17
    using ctype = typename HostGrid::ctype;
Praetorius, Simon's avatar
Praetorius, Simon committed
18 19 20

  public:
    // The codimension of this entitypointer wrt the host grid
21 22
    enum { codimension = HostGrid::dimension - mydim };
    enum { dimensionworld = HostGrid::dimensionworld };
Praetorius, Simon's avatar
Praetorius, Simon committed
23 24

    // select appropriate hostgrid geometry via typeswitch
25
    using HostLocalGeometry = typename HostGrid::Traits::template Codim<codimension>::LocalGeometry;
Praetorius, Simon's avatar
Praetorius, Simon committed
26 27

    /// entity in the host grid
28
    using HostGridEntity = typename HostGrid::Traits::template Codim<codimension>::Entity;
Praetorius, Simon's avatar
Praetorius, Simon committed
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148

    //! type of jacobian transposed
    using JacobianInverseTransposed = typename HostLocalGeometry::JacobianInverseTransposed;
    using JacobianTransposed = typename HostLocalGeometry::JacobianTransposed;


    /// Constructor from host entities
    MultiMeshLocalGeometry (const HostGridEntity& sourceEntity, const HostGridEntity& targetEntity)
    {
      hostLocalGeometries_.reserve(sourceEntity.level() - targetEntity.level());
      auto entity = sourceEntity;
      for (int l = sourceEntity.level(); l > targetEntity.level(); --l) {
        assert( entity.hasFather() );
        hostLocalGeometries_.push_back(entity.geometryInFather());
        if (l > targetEntity.level()+1)
          entity = entity.father();
      }
    }


    /// Return the element type identifier
    GeometryType type () const
    {
      return hostLocalGeometries_[0].type();
    }

    /// Return wether we have an affine mapping
    bool affine () const
    {
      return std::all_of(hostLocalGeometries_.begin(), hostLocalGeometries_.end(),
        [](auto const& localGeo) { return localGeo.affine(); });
    }

    /// Returns true if the point is in the current element
    bool checkInside (const FieldVector<ctype, mydim>& local) const
    {
      return hostLocalGeometries_[0].checkInside(local);
    }

    /// Return the number of corners of this element. Corners are numbered 0...n-1
    int corners () const
    {
      return hostLocalGeometries_[0].corners();
    }

    /// Access to coordinates of corners. Index is the number of the corner
    const FieldVector<ctype, coorddim> corner (int i) const
    {
      FieldVector<ctype, coorddim> c = hostLocalGeometries_[0].corner(i);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[j].global(c);
      return c;
    }


    /// Maps a local coordinate within reference element to global coordinate in element
    FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const
    {
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[j].global(c);
      return c;
    }

    /// Maps a global coordinate within the element to a local coordinate in its reference element
    FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const
    {
      auto c = hostLocalGeometries_.back().local(global);
      std::size_t n = hostLocalGeometries_.size()-1;
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[n-j].local(c);
      return c;
    }

    /// Return the factor appearing in the integral transformation formula
    ctype integrationElement (const FieldVector<ctype, mydim>& local) const
    {
      ctype dx = 1;
      auto c = local;
      for (std::size_t j = 0; j < hostLocalGeometries_.size(); ++j) {
        dx *= hostLocalGeometries_[j].integrationElement(c);
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return dx;
    }

    /// The Jacobian matrix of the mapping from the reference element to this element
    JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
    {
      auto J = hostLocalGeometries_[0].jacobianInverseTransposed(local);
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
        J.leftmultiply(hostLocalGeometries_[j].jacobianInverseTransposed(c));
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return J;
    }

    /// Return the transposed of the Jacobian
    JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const
    {
      auto J = hostLocalGeometries_[0].jacobianTransposed(local);
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
        J.rightmultiply(hostLocalGeometries_[j].jacobianTransposed(c));
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return J;
    }

  private:
    std::vector<HostLocalGeometry> hostLocalGeometries_;
  };

}  // end namespace Dune

#endif // DUNE_MULTIMESH_LOCALGEOMETRY_HH