geometry.hh

/*~--------------------------------------------------------------------------~*
 *  @@@@@@@@  @@           @@@@@@   @@@@@@@@ @@
 * /@@/////  /@@          @@////@@ @@////// /@@
 * /@@       /@@  @@@@@  @@    // /@@       /@@
 * /@@@@@@@  /@@ @@///@@/@@       /@@@@@@@@@/@@
 * /@@////   /@@/@@@@@@@/@@       ////////@@/@@
 * /@@       /@@/@@//// //@@    @@       /@@/@@
 * /@@       @@@//@@@@@@ //@@@@@@  @@@@@@@@ /@@
 * //       ///  //////   //////  ////////  //
 *
 * Copyright (c) 2016 Los Alamos National Laboratory, LLC
 * All rights reserved
 *~--------------------------------------------------------------------------~*/

#pragma once

#if !defined(__FLECSI_PRIVATE__)
#error Do not include this file directly!
#endif

#include "flecsi/util/geometry/point.hh"

#include <array>
#include <bitset>
#include <cassert>
#include <cmath>
#include <float.h>
#include <iostream>
#include <math.h>
#include <vector>

namespace flecsi {
namespace topo {

/*-----------------------------------------------------------------------------*
 * class ntree_geometry
 *-----------------------------------------------------------------------------*/
template<typename T, size_t D>
struct ntree_geometry {};

/*-----------------------------------------------------------------------------*
 * class ntree_geometry 1D specification
 *-----------------------------------------------------------------------------*/
template<typename T>
struct ntree_geometry<T, 1> {

  using point_t = util::point<T, 1>;
  using element_t = T;
  static constexpr element_t tol =
    std::numeric_limits<element_t>::epsilon() * 10.;

  static bool
  within(const point_t & origin, const point_t & center, element_t r1) {
    return util::distance(origin, center) - r1 <= tol;
  }

  static bool within_square(const point_t & origin,
    const point_t & center,
    element_t r1,
    element_t r2) {
    return util::distance(origin, center) <= std::max(r1, r2);
  }

  static bool
  within_box(const point_t & min, const point_t & max, const point_t & origin) {
    return origin[0] <= max[0] && origin[0] >= min[0];
  }

  static bool intersects_box_box(const point_t & min_b1,
    const point_t & max_b1,
    const point_t & min_b2,
    const point_t & max_b2) {
    return !((max_b1[0] < min_b2[0]) || (max_b2[0] < min_b1[0]));
  }

  static bool intersects_sphere_sphere(const point_t & c1,
    const element_t r1,
    const point_t & c2,
    const element_t r2) {
    return util::distance(c1, c2) - (r1 + r2) <= tol;
  }

  static bool intersects_sphere_box(const point_t & min,
    const point_t & max,
    const point_t & c,
    const element_t r) {
    point_t x = point_t(std::max(min[0], std::min(c[0], max[0])));
    element_t dist = util::distance(x, c);
    return dist - r <= tol;
  }

  bool box_MAC(const point_t & position_source,
    const point_t & position_sink,
    const point_t & box_source_min,
    const point_t & box_source_max,
    double macangle) {
    double dmax = util::distance(box_source_min, box_source_max);
    double disttoc = util::distance(position_sink, position_source);
    return dmax / disttoc - macangle <= tol;
  }
}; // class ntree_geometry specification for 1D

/*-----------------------------------------------------------------------------*
 * class ntree_geometry 2D specification
 *-----------------------------------------------------------------------------*/
template<typename T>
struct ntree_geometry<T, 2> {
  using point_t = util::point<T, 2>;
  using element_t = T;

  static constexpr element_t tol =
    std::numeric_limits<element_t>::epsilon() * 10.;

  static bool
  within(const point_t & origin, const point_t & center, element_t r1) {
    return util::distance(origin, center) - r1 <= tol;
  }

  static bool within_square(const point_t & origin,
    const point_t & center,
    element_t r1,
    element_t r2) {
    return util::distance(origin, center) <= std::max(r1, r2);
  }

  static bool
  within_box(const point_t & min, const point_t & max, const point_t & origin) {
    return origin[0] <= max[0] && origin[0] > min[0] && origin[1] <= max[1] &&
           origin[1] > min[1];
  }

  static bool intersects_box_box(const point_t & min_b1,
    const point_t & max_b1,
    const point_t & min_b2,
    const point_t & max_b2) {
    return !((max_b1[0] < min_b2[0] || max_b1[1] < min_b2[1]) ||
             (max_b2[0] < min_b1[0] || max_b2[1] < min_b1[1]));
  }

  static bool intersects_sphere_sphere(const point_t & c1,
    const element_t r1,
    const point_t & c2,
    const element_t r2) {
    return util::distance(c1, c2) - (r1 + r2) <= tol;
  }

  static bool intersects_sphere_box(const point_t & min,
    const point_t & max,
    const point_t & c,
    const element_t r) {
    point_t x = point_t(std::max(min[0], std::min(c[0], max[0])),
      std::max(min[1], std::min(c[1], max[1])));
    element_t dist = util::distance(x, c);
    return dist - r <= tol;
  }

  static bool box_MAC(const point_t & position_source,
    const point_t & position_sink,
    const point_t & box_source_min,
    const point_t & box_source_max,
    double macangle) {
    double dmax = util::distance(box_source_min, box_source_max);
    double disttoc = util::distance(position_sink, position_source);
    return dmax / disttoc < macangle;
  }
}; // class tree_geometry specification for 2D

/*-----------------------------------------------------------------------------*
 * class tree_geometry 3D specification
 *-----------------------------------------------------------------------------*/
template<typename T>
struct ntree_geometry<T, 3> {
  using point_t = util::point<T, 3>;
  using element_t = T;
  static constexpr element_t tol =
    std::numeric_limits<element_t>::epsilon() * 10.;

  static bool
  within(const point_t & origin, const point_t & center, element_t r1) {
    return util::distance(origin, center) - r1 <= tol;
  }

  inline static bool within_square(const point_t & origin,
    const point_t & center,
    const element_t & r1,
    const element_t & r2) {
    return util::distance(origin, center) <= std::max(r1, r2);
  }

  static bool
  within_box(const point_t & min, const point_t & max, const point_t & origin) {
    return origin[0] <= max[0] && origin[0] > min[0] && origin[1] <= max[1] &&
           origin[1] > min[1] && origin[2] <= max[2] && origin[2] > min[2];
  }

  inline static bool intersects_box_box(const point_t & min_b1,
    const point_t & max_b1,
    const point_t & min_b2,
    const point_t & max_b2) {
    return !((max_b1[0] < min_b2[0] || max_b1[1] < min_b2[1] ||
               max_b1[2] < min_b2[2]) ||
             (max_b2[0] < min_b1[0] || max_b2[1] < min_b1[1] ||
               max_b2[2] < min_b1[2]));
  }

  static bool intersects_sphere_sphere(const point_t & c1,
    const element_t r1,
    const point_t & c2,
    const element_t r2) {
    return (c2[0] - c1[0]) * (c2[0] - c1[0]) +
             (c2[1] - c1[1]) * (c2[1] - c1[1]) +
             (c2[2] - c1[2]) * (c2[2] - c1[2]) <=
           (r1 + r2) * (r1 + r2);
  }

  static bool intersects_sphere_box(const point_t & min,
    const point_t & max,
    const point_t & c,
    const element_t & r) {
    point_t x = point_t(c[0] < max[0] ? c[0] : max[0],
      c[1] < max[1] ? c[1] : max[1],
      c[2] < max[2] ? c[2] : max[2]);
    x = {x[0] < min[0] ? min[0] : x[0],
      x[1] < min[1] ? min[1] : x[1],
      x[2] < min[2] ? min[2] : x[2]};
    element_t dist = (x[0] - c[0]) * (x[0] - c[0]) +
                     (x[1] - c[1]) * (x[1] - c[1]) +
                     (x[2] - c[2]) * (x[2] - c[2]);
    return dist <= r * r;
  }

  static bool box_MAC(const point_t & position_source,
    const point_t & position_sink,
    const point_t & box_source_min,
    const point_t & box_source_max,
    double macangle) {
    double dmax = util::distance(box_source_min, box_source_max);
    double disttoc = util::distance(position_sink, position_source);
    return dmax / disttoc < macangle;
  }
}; // class ntree_geometry specification for 3D

} // namespace topo
} // namespace flecsi

/*~-------------------------------------------------------------------------~-*
 * Formatting options for vim.
 * vim: set tabstop=2 shiftwidth=2 expandtab :
 *~-------------------------------------------------------------------------~-*/