geometry.hpp

Go to the documentation of this file.

#pragma once
 
#include <Eigen/Core>
#include <array>
#include "math/differences.hpp"
#include "utils/assertion.hpp"
 
namespace precice {
namespace math {
 
namespace geometry {
 
enum ResultConstants {
  NO_INTERSECTION,
  INTERSECTION,
  TOUCHING,
  NOT_CONTAINED,
  CONTAINED,
};
 
bool lineIntersection(
    const Eigen::Ref<const Eigen::Vector2d> &a,
    const Eigen::Ref<const Eigen::Vector2d> &b,
    const Eigen::Ref<const Eigen::Vector2d> &c,
    const Eigen::Ref<const Eigen::Vector2d> &d,
    Eigen::Ref<Eigen::Vector2d> &            intersectionPoint);
 
ResultConstants segmentPlaneIntersection(
    const Eigen::Vector3d &pointOnPlane,
    const Eigen::Vector3d &planeNormal,
    const Eigen::Vector3d &firstPointSegment,
    const Eigen::Vector3d &secondPointSegment,
    Eigen::Vector3d &      intersectionPoint);
 
template <typename VECTORA_T, typename VECTORB_T, typename VECTORC_T>
bool between(
    const VECTORA_T &a,
    const VECTORB_T &b,
    const VECTORC_T &toCheck);
 
template <typename Derived>
bool collinear(
    const Eigen::MatrixBase<Derived> &a,
    const Eigen::MatrixBase<Derived> &b,
    const Eigen::MatrixBase<Derived> &c);
 
template <typename Derived>
static bool parallel(
    const Eigen::MatrixBase<Derived> &a,
    const Eigen::MatrixBase<Derived> &b,
    const Eigen::MatrixBase<Derived> &c,
    const Eigen::MatrixBase<Derived> &d);
 
double triangleArea(
    const Eigen::VectorXd &a,
    const Eigen::VectorXd &b,
    const Eigen::VectorXd &c);
 
double tetraVolume(
    const Eigen::Vector3d &a,
    const Eigen::Vector3d &b,
    const Eigen::Vector3d &c,
    const Eigen::Vector3d &d);
 
Eigen::Vector2d projectVector(
    const Eigen::Vector3d &vector,
    const int              indexDimensionToRemove);
 
template <class Derived>
int containedInHyperrectangle(
    const Eigen::MatrixBase<Derived> &sidelengths,
    const Eigen::MatrixBase<Derived> &center,
    const Eigen::MatrixBase<Derived> &testPoint);
 
// --------------------------------------------------------- HEADER DEFINITIONS
 
template <typename VECTORA_T, typename VECTORB_T, typename VECTORC_T>
bool between(
    const VECTORA_T &a,
    const VECTORB_T &b,
    const VECTORC_T &toCheck)
{
  if (not collinear(a, b, toCheck)) {
    return false;
  }
 
  if (a(0) != b(0)) {
    return (math::greaterEquals(toCheck(0), a(0)) && math::greaterEquals(b(0), toCheck(0))) ||
           (math::greaterEquals(a(0), toCheck(0)) && math::greaterEquals(toCheck(0), b(0)));
  } else {
    return (math::greaterEquals(toCheck(1), a(1)) && math::greaterEquals(b(1), toCheck(1))) ||
           (math::greaterEquals(a(1), toCheck(1)) && math::greaterEquals(toCheck(1), b(1)));
  }
}
 
template <typename Derived>
bool collinear(
    const Eigen::MatrixBase<Derived> &a,
    const Eigen::MatrixBase<Derived> &b,
    const Eigen::MatrixBase<Derived> &c)
{
  PRECICE_ASSERT(a.size() == b.size(), a.size(), b.size());
  PRECICE_ASSERT(a.size() == c.size(), a.size(), c.size());
  double triangleOutline = (b - a).norm() + (c - b).norm() + (a - c).norm();
  if (math::equals(triangleArea(a, b, c) / triangleOutline, 0.0)) {
    return true;
  }
  return false;
}
 
template <class Derived>
bool parallel(
    const Eigen::MatrixBase<Derived> &a,
    const Eigen::MatrixBase<Derived> &b,
    const Eigen::MatrixBase<Derived> &c,
    const Eigen::MatrixBase<Derived> &d)
{
  if (math::equals(triangleArea(a, b, c), 0.0) and math::equals(triangleArea(a, b, d), 0.0))
    return true;
 
  return false;
}
 
template <class Derived>
int containedInHyperrectangle(
    const Eigen::MatrixBase<Derived> &sidelengths,
    const Eigen::MatrixBase<Derived> &center,
    const Eigen::MatrixBase<Derived> &testPoint)
{
  int             dim      = sidelengths.size();
  Eigen::VectorXd toCenter = testPoint - center;
  toCenter                 = toCenter.cwiseAbs();
 
  double diff     = 0.0;
  bool   touching = false;
  for (int i = 0; i < dim; i++) {
    diff = 0.5 * sidelengths(i) - toCenter(i);
    if (math::greater(0.0, diff)) {
      return NOT_CONTAINED;
    }
    if (math::equals(diff, 0.0)) {
      touching = true;
    }
  }
  if (touching) {
    return TOUCHING;
  }
  return CONTAINED;
}
 
struct ConvexityResult {
  bool               convex;
  std::array<int, 4> vertexOrder;
};
 
ConvexityResult isConvexQuad(std::array<Eigen::VectorXd, 4> coords);
 
} // namespace geometry
} // namespace math
} // namespace precice