geometry.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_MMESH_GRID_GEOMETRY_HH
#define DUNE_MMESH_GRID_GEOMETRY_HH
 
// Dune includes
#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
#include <dune/geometry/affinegeometry.hh>
#include <dune/grid/common/geometry.hh>
 
// local includes
#include <dune/mmesh/grid/pointfieldvector.hh>
 
namespace Dune {
 
template <int mydim, int coorddim, class GridImp>
class MMeshGeometry {};
 
 
template <int md, class GridImp>
class MMeshGeometry<md, 2, GridImp>
    : public AffineGeometry<typename GridImp::ctype, md, 2> {
 public:
  static constexpr int mydim = md;
  static constexpr int coorddim = 2;
  typedef AffineGeometry<typename GridImp::ctype, mydim, coorddim> BaseType;
  typedef FieldVector<typename GridImp::ctype, coorddim> FVector;
  typedef typename GridImp::LeafIntersection Intersection;
 
  enum { dimension = GridImp::dimension };
  enum { dimensionworld = GridImp::dimensionworld };
  enum { coorddimension = coorddim };
  enum { mydimension = mydim };
 
  explicit MMeshGeometry()
      : BaseType(
            GeometryTypes::simplex(mydim),
            std::array<FVector, 3>({FVector({0.0, 0.0}), FVector({1.0, 0.0}),
                                    FVector({0.0, 1.0})})) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<0>& hostEntity)
      : BaseType(GeometryTypes::simplex(mydim), getVertices(hostEntity)) {}
 
  MMeshGeometry(const std::array<FVector, 3> points)
      : BaseType(GeometryTypes::simplex(mydim), points) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<1>& hostEntity)
      : BaseType(GeometryTypes::simplex(mydim), getVertices(hostEntity)) {}
 
  MMeshGeometry(int idx)
      : BaseType(GeometryTypes::simplex(mydim), getLocalVertices_(idx)) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<2>& hostEntity)
      : BaseType(
            GeometryTypes::simplex(mydim),
            std::array<FVector, 1>({makeFieldVector(hostEntity->point())})) {}
 
  const FVector circumcenter() const { return computeCircumcenter(*this); }
 
 private:
  static inline std::array<FVector, 3> getVertices(
      const typename GridImp::template HostGridEntity<0> hostEntity) {
    const auto& cgalIdx = hostEntity->info().cgalIndex;
 
    std::array<FVector, 3> vertices(
        {makeFieldVector(hostEntity->vertex(cgalIdx[0])->point()),
         makeFieldVector(hostEntity->vertex(cgalIdx[1])->point()),
         makeFieldVector(hostEntity->vertex(cgalIdx[2])->point())});
 
    return vertices;
  }
 
  static inline std::array<FVector, 2> getVertices(
      const typename GridImp::template HostGridEntity<1>& hostEntity) {
    const auto& cgalIdx = hostEntity.first->info().cgalIndex;
 
    auto facetIdx = MMeshImpl::cgalFacetToDuneFacet<
        2, typename GridImp::template HostGridEntity<1>>(hostEntity);
 
    std::array<FVector, 2> vertices;
    for (int k = 0; k < 2; ++k)
      vertices[k] = makeFieldVector(
          hostEntity.first
              ->vertex(
                  cgalIdx[MMeshImpl::ref<2>().subEntity(facetIdx, 1, k, 2)])
              ->point());
 
    return vertices;
  }
 
  static inline std::array<FVector, 2> getLocalVertices_(int k) {
    static const std::array<FVector, 3> local = {
        FVector({0.0, 0.0}), FVector({1.0, 0.0}), FVector({0.0, 1.0})};
 
    return {local[k == 2 ? 1 : 0], local[k == 0 ? 1 : 2]};
  }
};
 
 
template <int md, class GridImp>
class MMeshGeometry<md, 3, GridImp>
    : public AffineGeometry<typename GridImp::ctype, md, 3> {
 public:
  static constexpr int mydim = md;
  static constexpr int coorddim = 3;
  typedef AffineGeometry<typename GridImp::ctype, mydim, coorddim> BaseType;
  typedef typename GridImp::ctype ctype;
  typedef FieldVector<ctype, coorddim> FVector;
  typedef typename GridImp::LeafIntersection Intersection;
 
  enum { dimension = GridImp::dimension };
  enum { dimensionworld = GridImp::dimensionworld };
  enum { coorddimension = coorddim };
  enum { mydimension = mydim };
 
  explicit MMeshGeometry()
      : BaseType(GeometryTypes::simplex(mydim),
                 std::array<FVector, 4>(
                     {FVector({0.0, 0.0, 0.0}), FVector({1.0, 0.0, 0.0}),
                      FVector({0.0, 1.0, 0.0}), FVector({0.0, 0.0, 1.0})})) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<0>& hostEntity)
      : BaseType(GeometryTypes::simplex(mydim), getVertices<0>(hostEntity)) {}
 
  MMeshGeometry(const std::array<FVector, 4> points)
      : BaseType(GeometryTypes::simplex(mydim), points) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<1>& hostEntity)
      : BaseType(GeometryTypes::simplex(mydim), getVertices<1>(hostEntity)) {}
 
  MMeshGeometry(int idx)
      : BaseType(GeometryTypes::simplex(mydim), getLocalVertices_(idx)) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<2>& hostEntity)
      : BaseType(GeometryTypes::simplex(mydim), getVertices<2>(hostEntity)) {}
 
  MMeshGeometry(const typename GridImp::template HostGridEntity<3>& hostEntity)
      : BaseType(
            GeometryTypes::simplex(mydim),
            std::array<FVector, 1>({makeFieldVector(hostEntity->point())})) {}
 
  const FVector circumcenter() const { return computeCircumcenter(*this); }
 
 private:
  template <int codim, typename Enable = std::enable_if_t<codim == 0>>
  static inline std::array<FVector, 4> getVertices(
      const typename GridImp::template HostGridEntity<0>& hostEntity) {
    const auto& cgalIdx = hostEntity->info().cgalIndex;
 
    std::array<FVector, 4> vertices;
    for (int i = 0; i < 4; ++i)
      vertices[i] = makeFieldVector(hostEntity->vertex(cgalIdx[i])->point());
    return vertices;
  }
 
  template <int codim, typename Enable = std::enable_if_t<codim == 1>>
  static inline std::array<FVector, 3> getVertices(
      const typename GridImp::template HostGridEntity<1>& hostEntity) {
    const auto& cgalIdx = hostEntity.first->info().cgalIndex;
    const auto& cell = hostEntity.first;
 
    auto facetIdx = MMeshImpl::cgalFacetToDuneFacet<
        3, typename GridImp::template HostGridEntity<1>>(hostEntity);
 
    std::array<FVector, 3> vertices;
    for (int i = 0; i < 3; ++i)
      vertices[i] = makeFieldVector(
          cell->vertex(
                  cgalIdx[MMeshImpl::ref<3>().subEntity(facetIdx, 1, i, 3)])
              ->point());
    return vertices;
  }
 
  template <int codim, typename Enable = std::enable_if_t<codim == 2>>
  static inline std::array<FVector, 2> getVertices(
      const typename GridImp::template HostGridEntity<2>& hostEntity) {
    const auto& cgalIdx = hostEntity.first->info().cgalIndex;
    const auto& cell = hostEntity.first;
 
    auto edgeIdx = MMeshImpl::cgalEdgeToDuneEdge<
        3, typename GridImp::template HostGridEntity<2>>(hostEntity);
 
    std::array<FVector, 2> vertices;
    vertices[0] = makeFieldVector(
        cell->vertex(cgalIdx[MMeshImpl::ref<3>().subEntity(edgeIdx, 2, 0, 3)])
            ->point());
    vertices[1] = makeFieldVector(
        cell->vertex(cgalIdx[MMeshImpl::ref<3>().subEntity(edgeIdx, 2, 1, 3)])
            ->point());
 
    return vertices;
  }
 
  static inline std::array<FVector, 3> getLocalVertices_(int k) {
    static const std::array<FVector, 4> local = {
        FVector({0.0, 0.0, 0.0}), FVector({1.0, 0.0, 0.0}),
        FVector({0.0, 1.0, 0.0}), FVector({0.0, 0.0, 1.0})};
 
    return {local[k <= 2 ? 0 : 1], local[k <= 1 ? 1 : 2],
            local[k == 0 ? 2 : 3]};
  }
};
 
}  // namespace Dune
 
#endif