pointfieldvector.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_GRID_MMESHPOINTFIELDVECTOR_HH
#define DUNE_GRID_MMESHPOINTFIELDVECTOR_HH
 
// Dune includes
#include <dune/common/fvector.hh>
 
namespace Dune {
 
template <class Kernel>
static inline FieldVector<typename Kernel::RT, 2> makeFieldVector(
    const CGAL::Point_2<Kernel>& p) {
  return FieldVector<typename Kernel::RT, 2>({p.x(), p.y()});
}
 
template <class Kernel>
static inline FieldVector<typename Kernel::RT, 2> makeFieldVector(
    const CGAL::Vector_2<Kernel>& p) {
  return FieldVector<typename Kernel::RT, 2>({p.x(), p.y()});
}
 
template <class Kernel>
static inline FieldVector<typename Kernel::RT, 3> makeFieldVector(
    const CGAL::Point_3<Kernel>& p) {
  return FieldVector<typename Kernel::RT, 3>({p.x(), p.y(), p.z()});
}
 
template <class Kernel>
static inline FieldVector<typename Kernel::RT, 3> makeFieldVector(
    const CGAL::Vector_3<Kernel>& p) {
  return FieldVector<typename Kernel::RT, 3>({p.x(), p.y(), p.z()});
}
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel PointKernel;
 
template <typename ctype>
static inline auto makePoint(const Dune::FieldVector<ctype, 2>& v) {
  return CGAL::Point_2<PointKernel>(v[0], v[1]);
}
 
template <typename ctype>
static inline auto makePoint(const Dune::FieldVector<ctype, 3>& v) {
  return CGAL::Point_3<PointKernel>(v[0], v[1], v[2]);
}
 
template <class Geometry>
auto computeCircumcenter(const Geometry& geo) {
  static constexpr int mydim = Geometry::mydim;
 
  if constexpr (mydim == 3) {
    // obtain circumcenter by CGAL
    return makeFieldVector(
        CGAL::circumcenter(makePoint(geo.corner(0)), makePoint(geo.corner(1)),
                           makePoint(geo.corner(2)), makePoint(geo.corner(3))));
  }
 
  if constexpr (mydim == 2) {
    // obtain circumcenter by CGAL
    return makeFieldVector(CGAL::circumcenter(makePoint(geo.corner(0)),
                                              makePoint(geo.corner(1)),
                                              makePoint(geo.corner(2))));
  }
 
  if constexpr (mydim == 1) return 0.5 * (geo.corner(0) + geo.corner(1));
 
  if constexpr (mydim == 0) return geo.corner(0);
}
 
}  // namespace Dune
 
#endif