geomObjects.h

Go to the documentation of this file.

/*
* -------------------------------------------
* Copyright (c) 2021 - 2026 Prashant K. Jha
* -------------------------------------------
* PeriDEM https://github.com/prashjha/PeriDEM
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE)
*/
 
#ifndef GEOM_GEOMOBJECTS_H
#define GEOM_GEOMOBJECTS_H
 
#include <iostream>
#include <utility>
#include <memory>
#include <map>
 
#include "util/function.h"
#include "util/io.h"
#include "util/point.h" // definition of Point
#include "util/transformationFunctions.h"
 
namespace geom {
  
  const std::map<std::string, size_t> geom_type_to_dim = {
    {"circle", 2},
    {"ellipse", 2},
    {"square", 2},
    {"rectangle", 2},
    {"hexagon", 2},
    {"triangle", 2},
    {"drum2d", 2},
    {"sphere", 3},
    {"ellipsoid", 3},
    {"cube", 3},
    {"cuboid", 3},
    {"cylinder", 3},
    {"open_rect_channel_2d", 2},
    {"open_cuboid_channel_3d", 3},
    {"plane", 3},
    {"line", 1}
  };
 
  inline size_t getGeomTypeToDim(std::string type) {
    return geom_type_to_dim.at(type);
  };
 
  inline util::Point mapSimilarity(const util::Point &x, const util::Point &pivotDefault,
      const util::Point &t, double scale, double angle, const util::Point &axis,
      const util::Point *rotationPoint) {
    const util::Point &p = (rotationPoint != nullptr) ? *rotationPoint : pivotDefault;
    return p + util::rotate(x - p, angle, axis) * scale + t;
  }
 
struct BoxPartition {
 
util::Point d_xc; // centroid of box
std::pair<util::Point, util::Point> d_box; // two corner points
double d_r; // radius of circle inscribing box
std::vector<size_t> d_nodes; // id of nodes in this box
 
BoxPartition() : d_xc(util::Point()),
                 d_box({util::Point(), util::Point()}),
                 d_r(0.) {}
 
bool isNear(const BoxPartition &box, const double &tol) const {
  auto dx = d_xc - box.d_xc;
  return dx.length() < d_r + box.d_r + tol;
}
 
bool isNear(const util::Point &x, const double &tol) const {
 
  if (util::isLess(x.d_x, d_box.first.d_x - tol) or
      util::isLess(x.d_y, d_box.first.d_y - tol) or
      util::isLess(x.d_z, d_box.first.d_z - tol) or
      util::isGreater(x.d_x, d_box.second.d_x + tol) or
      util::isGreater(x.d_y, d_box.second.d_y + tol) or
      util::isGreater(x.d_z, d_box.second.d_z + tol))
    return false;
 
  return true;
}
 
void addNode(const size_t &i) {
  for (const auto &j: d_nodes)
    if (j == i)
      return;
 
  d_nodes.push_back(i);
}
};
 
 
  struct WallContactHit {
    bool active = false;
    double signed_gap = 0.;
    util::Point outward_n;
    util::Point closest;
  };
 
  class GeomObject {
 
public:
    const std::string d_name;
 
    const std::string d_description;
 
    std::vector<std::string> d_tags;
 
  public:
    explicit GeomObject(std::string name = "",
                        std::string description = "")
      : d_name(name), d_description(description) {
    };
 
    virtual double volume() const { return 0.; };
 
    virtual util::Point center() const { return {}; };
 
    virtual std::pair<util::Point, util::Point> box() const {
      return {util::Point(), util::Point()};
    };
 
    virtual std::pair<util::Point, util::Point>
    box(const double &tol) const {
      return {util::Point(), util::Point()};
    };
 
    virtual double inscribedRadius() const { return 0.; };
 
    virtual double boundingRadius() const { return 0.; };
 
    virtual bool isInside(const util::Point &x) const { return false; };
 
    virtual bool
    isOutside(const util::Point &x) const { return false; };
 
    virtual bool isNear(const util::Point &x, const double &tol) const {
      return false;
    };
 
    virtual bool isNearBoundary(const util::Point &x, const double &tol,
                                const bool &within) const {
      return false;
    };
 
    virtual bool
    doesIntersect(const util::Point &x) const { return false; };
 
    virtual bool
    isInside(const std::pair<util::Point, util::Point> &box) const {
      return false;
    };
 
    virtual bool
    isOutside(const std::pair<util::Point, util::Point> &box) const {
      return false;
    };
 
    virtual bool isNear(const std::pair<util::Point, util::Point> &box,
                        const double &tol) const {
      return false;
    };
 
    virtual bool
    doesIntersect(
      const std::pair<util::Point, util::Point> &box) const {
      return false;
    };
 
    virtual std::string printStr(int nt, int lvl) const {
      auto tabS = util::io::getTabS(nt);
      std::ostringstream oss;
 
      oss << tabS << "------- GeomObject --------" << std::endl
          << std::endl;
      oss << tabS << "Base geometrical object" << std::endl;
      oss << tabS << "Base implementation of GeomObject." << std::endl;
      oss << tabS << "Name of GeomObject = " << d_name << std::endl;
      oss << tabS << "Description of GeomObject = " << d_description
          << std::endl;
 
      return oss.str();
    };
 
    virtual void print(int nt, int lvl) const {
      std::cout << printStr(nt, lvl);
    };
 
    virtual void print() const { print(0, 0); };
 
    virtual void transform(const util::Point &translation, const double &scale, 
                         const double &angle, const util::Point &axis,
                         const util::Point *rotationPoint = nullptr) {
      // Base implementation does nothing
    }
 
    virtual bool wallContactQuery(const util::Point &x,
                                  WallContactHit &hit) const {
      (void)x;
      hit = WallContactHit();
      return false;
    }
  };
 
  class NullGeomObject : public GeomObject {
  public:
    NullGeomObject(std::string description = "") : GeomObject("null", description) {
    };
 
    NullGeomObject(const NullGeomObject &other) 
      : GeomObject(other.d_name, other.d_description) {
      d_tags = other.d_tags;
    }
 
    NullGeomObject& operator=(const NullGeomObject& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis,
                  const util::Point *rotationPoint) override {
      (void)translation;
      (void)scale;
      (void)angle;
      (void)axis;
      (void)rotationPoint;
    }
 
    virtual std::string printStr(int nt, int lvl) const override {
      auto tabS = util::io::getTabS(nt);
      std::ostringstream oss;
 
      oss << tabS << "------- NullGeomObject --------" << std::endl
          << std::endl;
      oss << tabS << "Name of GeomObject = " << d_name << std::endl;
 
      return oss.str();
    };
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Plane : public GeomObject {
  public:
    util::Point d_n; 
    util::Point d_p; 
 
    Plane()
      : GeomObject("plane", ""),
        d_n(util::Point(0., 1., 0.)),
        d_p(util::Point()) {}
 
    Plane(util::Point n, util::Point p, std::string description = "")
      : GeomObject("plane", std::move(description)), d_n(n), d_p(p) {
      const double ln = d_n.length();
      if (ln > 1.e-16)
        d_n = d_n / ln;
    }
 
    Plane(const Plane &other)
      : GeomObject(other.d_name, other.d_description),
        d_n(other.d_n),
        d_p(other.d_p) {
      d_tags = other.d_tags;
    }
 
    Plane &operator=(const Plane &other) {
      if (this != &other) {
        d_tags = other.d_tags;
        d_n = other.d_n;
        d_p = other.d_p;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale,
                   const double &angle, const util::Point &axis,
                   const util::Point *rotationPoint) override {
      (void)scale;
      const util::Point pivot =
          (rotationPoint != nullptr) ? *rotationPoint : d_p;
      d_p = mapSimilarity(d_p, pivot, translation, 1., angle, axis, rotationPoint);
      d_n = util::rotate(d_n, angle, axis);
      const double ln = d_n.length();
      if (ln > 1.e-16)
        d_n = d_n / ln;
    }
 
    double volume() const override { return 0.; }
    util::Point center() const override { return d_p; }
    std::pair<util::Point, util::Point> box() const override {
      return {d_p, d_p};
    }
    std::pair<util::Point, util::Point> box(const double &tol) const override {
      return {d_p - util::Point(tol, tol, tol), d_p + util::Point(tol, tol, tol)};
    }
    double inscribedRadius() const override { return 0.; }
    double boundingRadius() const override { return 0.; }
 
    bool isInside(const util::Point &x) const override {
      return (x - d_p) * d_n < 0.;
    }
    bool isOutside(const util::Point &x) const override {
      return !isInside(x);
    }
    bool isNear(const util::Point &x, const double &tol) const override {
      return std::abs((x - d_p) * d_n) <= tol;
    }
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool & /*within*/) const override {
      return isNear(x, tol);
    }
    bool doesIntersect(const util::Point &x) const override {
      return isNear(x, 1.e-12);
    }
 
    bool wallContactQuery(const util::Point &x,
                          WallContactHit &hit) const override;
 
    std::string printStr(int nt, int lvl) const override {
      auto tabS = util::io::getTabS(nt);
      std::ostringstream oss;
      oss << tabS << "------- Plane --------\n\n";
      oss << tabS << "Normal = " << d_n.printStr(0, lvl) << "\n";
      oss << tabS << "Point  = " << d_p.printStr(0, lvl) << "\n";
      return oss.str();
    }
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    }
    void print() const override { print(0, 0); }
  };
 
  class Line : public GeomObject {
  public: 
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_L;
 
    double d_r;
 
  public:
    Line()
      : GeomObject("line", ""),
        d_r(0.),
        d_L(0.),
        d_x(util::Point()),
        d_vertices({util::Point(), util::Point()}) {
    };
 
    Line(util::Point x1, util::Point x2, std::string description = "")
      : GeomObject("line", description),
        d_L((x1 - x2).length()),
        d_r(0.5 * (x1 - x2).length()),
        d_x(0.5 * (x1 + x2)),
        d_vertices({x1, x2}) {
    };
 
    Line(double L, util::Point x = util::Point(),
         std::string description = "")
      : GeomObject("line", description),
        d_L(L),
        d_r(0.5 * L),
        d_x(x),
        d_vertices({
          x + util::Point(-0.5 * L, 0., 0.),
          x + util::Point(0.5 * L, 0., 0.)
        }) {
    };
 
    Line(const Line &other)
      : GeomObject(other.d_name, other.d_description),
        d_r(other.d_r),
        d_L(other.d_L),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Line& operator=(const Line& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_r = other.d_r;
        d_L = other.d_L;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the length
      d_L *= scale;
      d_r *= scale;
 
      const util::Point c0 = d_x;
      // Rotate and translate vertices
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    bool wallContactQuery(const util::Point &x,
                          WallContactHit &hit) const override;
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Triangle : public GeomObject {
  public: 
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_r;
 
    util::Point d_a;
 
  public:
    Triangle()
      : GeomObject("triangle", ""),
        d_r(0.),
        d_a(util::Point(1., 0., 0.)),
        d_x(util::Point()),
        d_vertices(std::vector<util::Point>(3, util::Point())) {
    };
 
    Triangle(double r, util::Point x = util::Point(0., 0., 0.),
             util::Point a = util::Point(1., 0., 0.),
             std::string description = "uniform")
      : GeomObject("triangle", description),
        d_r(r),
        d_a(a),
        d_x(x),
        d_vertices(std::vector<util::Point>(3, util::Point())) {
      // generate vertices
      auto rotate_axis = util::Point(0., 0., 1.); // z-axis
      for (int i = 0; i < 3; i++) {
        d_vertices[i] = d_x + d_r * util::rotate(
                          d_a, i * 2. * M_PI / 3., rotate_axis);
      }
    };
 
    Triangle(util::Point v0, util::Point v1, util::Point v2,
             std::string description = "vertices")
      : GeomObject("triangle", description),
        d_r(0.),
        d_a(util::Point(1., 0., 0.)),
        d_x(util::Point()),
        d_vertices({v0, v1, v2}) {
      d_x = (v0 + v1 + v2) / 3.0;
      d_r = (v0 - d_x).length();
      {
        const double r1 = (v1 - d_x).length();
        const double r2 = (v2 - d_x).length();
        if (r1 > d_r)
          d_r = r1;
        if (r2 > d_r)
          d_r = r2;
      }
      d_a = v0 - d_x;
      if (d_a.length() > 0.)
        d_a = d_a / d_a.length();
    };
 
    Triangle(const Triangle &other)
      : GeomObject(other.d_name, other.d_description),
        d_r(other.d_r),
        d_a(other.d_a),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Triangle& operator=(const Triangle& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_r = other.d_r;
        d_a = other.d_a;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the radius
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
      d_a = util::rotate(d_a, angle, axis);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Square : public GeomObject {
 
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_L;
 
    double d_r;
 
  public:
    Square()
      : GeomObject("square", ""),
        d_L(0.),
        d_r(0.),
        d_x(util::Point()),
        d_vertices({
          util::Point(), util::Point(),
          util::Point(), util::Point()
        }) {
    };
 
    Square(double L, util::Point x = util::Point(),
           std::string description = "")
      : GeomObject("square", description),
        d_L(L),
        d_r(L / std::sqrt(2)),
        d_x(x),
        d_vertices({
          x + util::Point(-0.5 * L, -0.5 * L, 0.),
          x + util::Point(0.5 * L, -0.5 * L, 0.),
          x + util::Point(0.5 * L, 0.5 * L, 0.),
          x + util::Point(-0.5 * L, 0.5 * L, 0.)
        }) {
    };
 
    Square(util::Point x1, util::Point x2, std::string description = "")
      : GeomObject("square", description),
        d_L((x1 - x2).length() * (1. / std::sqrt(2))),
        d_r((x1 - x2).length() / 2.),
        d_x(0.5 * (x1 + x2)) {
      d_vertices.push_back(x1);
      d_vertices.push_back(util::Point(x1.d_x + d_L, x1.d_y, x1.d_z));
      d_vertices.push_back(x2);
      d_vertices.push_back(util::Point(x2.d_x - d_L, x2.d_y, x2.d_z));
    };
 
    Square(const Square &other)
      : GeomObject(other.d_name, other.d_description),
        d_L(other.d_L),
        d_r(other.d_r),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Square& operator=(const Square& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_L = other.d_L;
        d_r = other.d_r;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_L *= scale;
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Rectangle : public GeomObject {
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_Lx;
 
    double d_Ly;
 
    double d_r;
    
  public:
    Rectangle()
      : GeomObject("rectangle", ""),
        d_Lx(0.),
        d_Ly(0.),
        d_r(0.),
        d_x(util::Point()),
        d_vertices({
          util::Point(), util::Point(),
          util::Point(), util::Point()
        }) {
    };
 
    Rectangle(double Lx, double Ly, util::Point x = util::Point(),
              std::string description = "")
      : GeomObject("rectangle", description),
        d_Lx(Lx),
        d_Ly(Ly),
        d_r(0.5 * std::sqrt(std::pow(Lx, 2) + std::pow(Ly, 2))),
        d_x(x) {
      d_vertices.push_back(
        x + util::Point(-0.5 * d_Lx, -0.5 * d_Ly, 0.));
      d_vertices.push_back(
        x + util::Point(0.5 * d_Lx, -0.5 * d_Ly, 0.));
      d_vertices.push_back(x + util::Point(0.5 * d_Lx, 0.5 * d_Ly, 0.));
      d_vertices.push_back(
        x + util::Point(-0.5 * d_Lx, 0.5 * d_Ly, 0.));
    };
 
    Rectangle(const Rectangle &other)
      : GeomObject(other.d_name, other.d_description),
        d_Lx(other.d_Lx),
        d_Ly(other.d_Ly),
        d_r(other.d_r),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Rectangle(util::Point x1, util::Point x2,
              std::string description = "")
      : GeomObject("rectangle", description),
        d_Lx(x2.d_x - x1.d_x),
        d_Ly(x2.d_y - x1.d_y),
        d_r(0.5 *
            std::sqrt(std::pow(d_Lx, 2) + std::pow(d_Ly, 2))),
        d_x(0.5 * (x1 + x2)) {
      d_vertices.push_back(x1);
      d_vertices.push_back(util::Point(x1.d_x + d_Lx, x1.d_y, x1.d_z));
      d_vertices.push_back(x2);
      d_vertices.push_back(util::Point(x2.d_x - d_Lx, x2.d_y, x2.d_z));
    };
 
    Rectangle& operator=(const Rectangle& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_Lx = other.d_Lx;
        d_Ly = other.d_Ly;
        d_r = other.d_r;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
 
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_Lx *= scale;
      d_Ly *= scale;
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    bool wallContactQuery(const util::Point &x,
                          WallContactHit &hit) const override;
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Hexagon : public GeomObject {
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_r;
 
    util::Point d_a;
 
  public:
    Hexagon()
      : GeomObject("hexagon", ""),
        d_r(0.),
        d_a(util::Point(1., 0., 0.)),
        d_x(util::Point()),
        d_vertices(std::vector<util::Point>(6, util::Point())) {
    };
 
    Hexagon(double r, util::Point x = util::Point(0., 0., 0.),
            util::Point a = util::Point(1., 0., 0.),
            std::string description = "")
      : GeomObject("hexagon", ""),
        d_r(r),
        d_a(a),
        d_x(x),
        d_vertices(std::vector<util::Point>(6, util::Point())) {
      // generate vertices
      auto rotate_axis = util::Point(0., 0., 1.); // z-axis
      for (int i = 0; i < 6; i++) {
        d_vertices[i] = d_x + d_r * util::rotate(
                          d_a, i * M_PI / 3., rotate_axis);
      }
    };
 
    Hexagon(const Hexagon &other)
      : GeomObject(other.d_name, other.d_description),
        d_r(other.d_r),
        d_a(other.d_a),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Hexagon& operator=(const Hexagon& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_r = other.d_r;
        d_a = other.d_a;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the radius
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
      d_a = util::rotate(d_a, angle, axis);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Drum2D : public GeomObject {
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_w;
 
    double d_r;
 
    util::Point d_a;
 
  public:
    Drum2D()
      : GeomObject("drum2d", ""),
        d_r(0.),
        d_w(0.),
        d_a(util::Point(1., 0., 0.)),
        d_x(util::Point()),
        d_vertices(std::vector<util::Point>(6, util::Point())) {
    };
 
    Drum2D(double r, double w, util::Point x = util::Point(0., 0., 0.),
           util::Point a = util::Point(1., 0., 0.),
           std::string description = "")
      : GeomObject("drum2d", description),
        d_r(r),
        d_w(w),
        d_a(a),
        d_x(x),
        d_vertices(std::vector<util::Point>(6, util::Point())) {
      // generate vertices in the plane perpendicular to rotate_axis (z for axis a in xy)
      const util::Point rotate_axis(0., 0., 1.);
 
      d_vertices[0] = d_x + d_w * d_a;
      d_vertices[3] = d_x - d_w * d_a;
 
      d_vertices[1] = d_x + d_r * util::rotate(d_a, M_PI / 3., rotate_axis);
      /* Option 2 used (v1 - k*z_hat), which is not in the drum plane for a in xy and
       * breaks coplanarity (e.g. z != 0 on some corners) -> degenerate nodal volumes. */
      d_vertices[2] = d_x + d_r * util::rotate(d_a, 2. * M_PI / 3., rotate_axis);
 
      d_vertices[4] = d_x + d_r * util::rotate(-1. * d_a, M_PI / 3., rotate_axis);
      d_vertices[5] = d_x + d_r * util::rotate(-1. * d_a, 2. * M_PI / 3., rotate_axis);
    };
 
    Drum2D& operator=(const Drum2D& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_r = other.d_r;
        d_w = other.d_w;
        d_a = other.d_a;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_r *= scale;
      d_w *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
      d_a = util::rotate(d_a, angle, axis);
    }  
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Cube : public GeomObject {
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_L;
 
    double d_r;
 
  public:
    Cube()
      : GeomObject("cube", ""),
        d_L(0.),
        d_r(0.),
        d_x(util::Point()),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
    };
 
    Cube(double L, util::Point x = util::Point(),
         std::string description = "")
      : GeomObject("cube", description),
        d_L(L),
        d_r(0.5 * std::sqrt(3.) * L),
        d_x(x),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
      // numbering assuming x - left-right, y - back-front, z - bottom-top
      // counterclockwise in bottom plane (i.e.,
      // 0 - left-back-bottom, 1 - right-back-bottom,
      // 2 - right-front-bottom, 3 - left-front-bottom
      // similarly, in top plane
      d_vertices[0] =
          d_x + util::Point(-0.5 * d_L, -0.5 * d_L, -0.5 * d_L);
      d_vertices[1] =
          d_x + util::Point(0.5 * d_L, -0.5 * d_L, -0.5 * d_L);
      d_vertices[2] =
          d_x + util::Point(0.5 * d_L, 0.5 * d_L, -0.5 * d_L);
      d_vertices[3] =
          d_x + util::Point(-0.5 * d_L, 0.5 * d_L, -0.5 * d_L);
 
      d_vertices[4] =
          d_x + util::Point(-0.5 * d_L, -0.5 * d_L, 0.5 * d_L);
      d_vertices[5] =
          d_x + util::Point(0.5 * d_L, -0.5 * d_L, 0.5 * d_L);
      d_vertices[6] =
          d_x + util::Point(0.5 * d_L, 0.5 * d_L, 0.5 * d_L);
      d_vertices[7] =
          d_x + util::Point(-0.5 * d_L, 0.5 * d_L, 0.5 * d_L);
    };
 
    Cube(const Cube &other)
      : GeomObject(other.d_name, other.d_description),
        d_L(other.d_L),
        d_r(other.d_r),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Cube(util::Point x1, util::Point x2, std::string description = "")
      : GeomObject("cube", description),
        d_L((x2 - x1).length() / std::sqrt(3.)),
        d_r(0.5 * (x2 - x1).length()),
        d_x(0.5 * (x1 + x2)),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
      // numbering assuming x - left-right, y - back-front, z - bottom-top
      // counterclockwise in bottom plane (i.e.,
      // 0 - left-back-bottom, 1 - right-back-bottom,
      // 2 - right-front-bottom, 3 - left-front-bottom
      // similarly, in top plane
      d_vertices[0] =
          d_x + util::Point(-0.5 * d_L, -0.5 * d_L, -0.5 * d_L);
      d_vertices[1] =
          d_x + util::Point(0.5 * d_L, -0.5 * d_L, -0.5 * d_L);
      d_vertices[2] =
          d_x + util::Point(0.5 * d_L, 0.5 * d_L, -0.5 * d_L);
      d_vertices[3] =
          d_x + util::Point(-0.5 * d_L, 0.5 * d_L, -0.5 * d_L);
 
      d_vertices[4] =
          d_x + util::Point(-0.5 * d_L, -0.5 * d_L, 0.5 * d_L);
      d_vertices[5] =
          d_x + util::Point(0.5 * d_L, -0.5 * d_L, 0.5 * d_L);
      d_vertices[6] =
          d_x + util::Point(0.5 * d_L, 0.5 * d_L, 0.5 * d_L);
      d_vertices[7] =
          d_x + util::Point(-0.5 * d_L, 0.5 * d_L, 0.5 * d_L);
    };
 
    Cube& operator=(const Cube& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_L = other.d_L;
        d_r = other.d_r;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_L *= scale;
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Cuboid : public GeomObject {
  public:
    std::vector<util::Point> d_vertices;
 
    util::Point d_x;
 
    double d_Lx;
 
    double d_Ly;
 
    double d_Lz;
 
    double d_r;
 
  public:
    Cuboid()
      : GeomObject("cuboid", ""),
        d_Lx(0.),
        d_Ly(0.),
        d_Lz(0.),
        d_r(0.),
        d_x(util::Point()),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
    };
 
    Cuboid(double Lx, double Ly, double Lz,
           util::Point x = util::Point(),
           std::string description = "")
      : GeomObject("cuboid", description),
        d_Lx(Lx),
        d_Ly(Ly),
        d_Lz(Lz),
        d_r(std::sqrt(std::pow(Lx, 2) + std::pow(Ly, 2) +
                      std::pow(Lz, 2))),
        d_x(x),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
      // numbering assuming x - left-right, y - back-front, z - bottom-top
      // counterclockwise in bottom plane (i.e.,
      // 0 - left-back-bottom, 1 - right-back-bottom,
      // 2 - right-front-bottom, 3 - left-front-bottom
      // similarly, in top plane
      d_vertices[0] =
          d_x + util::Point(-0.5 * d_Lx, -0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[1] =
          d_x + util::Point(0.5 * d_Lx, -0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[2] =
          d_x + util::Point(0.5 * d_Lx, 0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[3] =
          d_x + util::Point(-0.5 * d_Lx, 0.5 * d_Ly, -0.5 * d_Lz);
 
      d_vertices[4] =
          d_x + util::Point(-0.5 * d_Lx, -0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[5] =
          d_x + util::Point(0.5 * d_Lx, -0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[6] =
          d_x + util::Point(0.5 * d_Lx, 0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[7] =
          d_x + util::Point(-0.5 * d_Lx, 0.5 * d_Ly, 0.5 * d_Lz);
    };
 
    Cuboid(const Cuboid &other)
      : GeomObject(other.d_name, other.d_description),
        d_Lx(other.d_Lx),
        d_Ly(other.d_Ly),
        d_Lz(other.d_Lz),
        d_r(other.d_r),
        d_x(other.d_x),
        d_vertices(other.d_vertices) {
      d_tags = other.d_tags;
    }
 
    Cuboid(util::Point x1, util::Point x2, std::string description = "")
      : GeomObject("cuboid", description),
        d_Lx(x2.d_x - x1.d_x),
        d_Ly(x2.d_y - x1.d_y),
        d_Lz(x2.d_z - x1.d_z),
        d_r(std::sqrt(
            std::pow(x2.d_x - x1.d_x, 2)
            + std::pow(x2.d_y - x1.d_y, 2)
            + std::pow(x2.d_z - x1.d_z, 2)
          )
        ),
        d_x(0.5 * (x1 + x2)),
        d_vertices(std::vector<util::Point>(8, util::Point())) {
      // numbering assuming x - left-right, y - back-front, z - bottom-top
      // counterclockwise in bottom plane (i.e.,
      // 0 - left-back-bottom, 1 - right-back-bottom,
      // 2 - right-front-bottom, 3 - left-front-bottom
      // similarly, in top plane
      d_vertices[0] =
          d_x + util::Point(-0.5 * d_Lx, -0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[1] =
          d_x + util::Point(0.5 * d_Lx, -0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[2] =
          d_x + util::Point(0.5 * d_Lx, 0.5 * d_Ly, -0.5 * d_Lz);
      d_vertices[3] =
          d_x + util::Point(-0.5 * d_Lx, 0.5 * d_Ly, -0.5 * d_Lz);
 
      d_vertices[4] =
          d_x + util::Point(-0.5 * d_Lx, -0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[5] =
          d_x + util::Point(0.5 * d_Lx, -0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[6] =
          d_x + util::Point(0.5 * d_Lx, 0.5 * d_Ly, 0.5 * d_Lz);
      d_vertices[7] =
          d_x + util::Point(-0.5 * d_Lx, 0.5 * d_Ly, 0.5 * d_Lz);
    };
 
    Cuboid& operator=(const Cuboid& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_Lx = other.d_Lx;
        d_Ly = other.d_Ly;
        d_Lz = other.d_Lz;
        d_r = other.d_r;
        d_x = other.d_x;
        d_vertices = other.d_vertices;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_Lx *= scale;
      d_Ly *= scale;
      d_Lz *= scale;
      d_r *= scale;
 
      const util::Point c0 = d_x;
      for (auto &v : d_vertices) {
        v = mapSimilarity(v, c0, translation, scale, angle, axis, rotationPoint);
      }
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Circle : public GeomObject {
  public:
    util::Point d_x;
 
    double d_r;
 
  public:
    Circle()
      : GeomObject("circle", ""),
        d_x(util::Point()),
        d_r(0.) {
    };
 
    Circle(double r, util::Point x = util::Point(),
           std::string description = "")
      : GeomObject("circle", description),
        d_x(x),
        d_r(r) {
    };
 
    Circle(const Circle &other)
      : GeomObject(other.d_name, other.d_description),
        d_x(other.d_x),
        d_r(other.d_r) {
      d_tags = other.d_tags;
    }
 
    Circle& operator=(const Circle& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_x = other.d_x;
        d_r = other.d_r;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      const util::Point c0 = d_x;
      d_r *= scale;
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Ellipse : public GeomObject {
  public:
    util::Point d_x;
 
    double d_a;
 
    double d_b;
 
    double d_theta;
 
  public:
    Ellipse()
      : GeomObject("ellipse", ""),
        d_x(util::Point()),
        d_a(0.),
        d_b(0.),
        d_theta(0.) {};
 
    Ellipse(double a, double b, double theta, util::Point x = util::Point(),
            std::string description = "")
      : GeomObject("ellipse", description),
        d_x(x),
        d_a(a),
        d_b(b),
        d_theta(theta) {};
 
    Ellipse(const Ellipse &other)
      : GeomObject(other.d_name, other.d_description),
        d_x(other.d_x),
        d_a(other.d_a),
        d_b(other.d_b),
        d_theta(other.d_theta) {
      d_tags = other.d_tags;
    }
 
    Ellipse &operator=(const Ellipse &other) {
      if (this != &other) {
        d_tags = other.d_tags;
        d_x = other.d_x;
        d_a = other.d_a;
        d_b = other.d_b;
        d_theta = other.d_theta;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, const double &angle,
                   const util::Point &axis, const util::Point *rotationPoint) override {
      const util::Point c0 = d_x;
      // Assumption: axis of rotation is the z-axis; d_theta is in-plane on top of world map
      d_a *= scale;
      d_b *= scale;
      d_theta += angle;
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
    util::Point center() const override;
    std::pair<util::Point, util::Point> box() const override;
    std::pair<util::Point, util::Point> box(const double &tol) const override;
    double inscribedRadius() const override;
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
    bool isOutside(const util::Point &x) const override;
    bool isNear(const util::Point &x, const double &tol) const override;
    bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override;
    bool doesIntersect(const util::Point &x) const override;
 
    bool isInside(const std::pair<util::Point, util::Point> &box) const override;
    bool isOutside(const std::pair<util::Point, util::Point> &box) const override;
    bool isNear(const std::pair<util::Point, util::Point> &box, const double &tol) const override;
    bool doesIntersect(const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
    void print(int nt, int lvl) const override { std::cout << printStr(nt, lvl); };
    void print() const override { print(0, 0); };
  };
 
  class Sphere : public GeomObject {
  public:
    util::Point d_x;
 
    double d_r;
 
  public:
    Sphere()
      : GeomObject("sphere", ""),
        d_x(util::Point()),
        d_r(0.) {
    };
 
    Sphere(double r, util::Point x = util::Point(),
           std::string description = "")
      : GeomObject("sphere", description),
        d_x(x),
        d_r(r) {
    };
 
    Sphere(const Sphere &other)
      : GeomObject(other.d_name, other.d_description),
        d_x(other.d_x),
        d_r(other.d_r) {
      d_tags = other.d_tags;
    }
 
    Sphere& operator=(const Sphere& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_x = other.d_x;
        d_r = other.d_r;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      const util::Point c0 = d_x;
      d_r *= scale;
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
 
  class Ellipsoid : public GeomObject {
  public:
    util::Point d_x;
    double d_a;
    double d_b;
    double d_c;
    util::Point d_axis;
    double d_theta;
 
  public:
    Ellipsoid()
      : GeomObject("ellipsoid", ""),
        d_x(util::Point()),
        d_a(0.),
        d_b(0.),
        d_c(0.),
        d_axis(0., 0., 1.),
        d_theta(0.) {};
 
    Ellipsoid(double a, double b, double c, util::Point x = util::Point(),
              std::string description = "")
      : GeomObject("ellipsoid", description),
        d_x(x),
        d_a(a),
        d_b(b),
        d_c(c),
        d_axis(0., 0., 1.),
        d_theta(0.) {};
 
    Ellipsoid(double a, double b, double c, double theta, util::Point x, util::Point axis,
              std::string description = "")
      : GeomObject("ellipsoid", description),
        d_x(x),
        d_a(a),
        d_b(b),
        d_c(c),
        d_axis(axis),
        d_theta(theta) {
      
      if (std::abs(theta) > 1.0e-12 and axis.length() > 1.0e-12) {
        d_axis = axis / axis.length();
      }
    }
 
    Ellipsoid(const Ellipsoid &other)
      : GeomObject(other.d_name, other.d_description),
        d_x(other.d_x),
        d_a(other.d_a),
        d_b(other.d_b),
        d_c(other.d_c),
        d_axis(other.d_axis),
        d_theta(other.d_theta) {
      d_tags = other.d_tags;
    }
 
    Ellipsoid &operator=(const Ellipsoid &other) {
      if (this != &other) {
        d_tags = other.d_tags;
        d_x = other.d_x;
        d_a = other.d_a;
        d_b = other.d_b;
        d_c = other.d_c;
        d_axis = other.d_axis;
        d_theta = other.d_theta;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, const double &angle,
                   const util::Point &axis, const util::Point *rotationPoint) override;
 
    double volume() const override;
    util::Point center() const override;
    std::pair<util::Point, util::Point> box() const override;
    std::pair<util::Point, util::Point> box(const double &tol) const override;
    double inscribedRadius() const override;
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
    bool isOutside(const util::Point &x) const override;
    bool isNear(const util::Point &x, const double &tol) const override;
    bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override;
    bool doesIntersect(const util::Point &x) const override;
 
    bool isInside(const std::pair<util::Point, util::Point> &box) const override;
    bool isOutside(const std::pair<util::Point, util::Point> &box) const override;
    bool isNear(const std::pair<util::Point, util::Point> &box, const double &tol) const override;
    bool doesIntersect(const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
    void print(int nt, int lvl) const override { std::cout << printStr(nt, lvl); };
    void print() const override { print(0, 0); };
  };
 
  class Cylinder : public GeomObject {
  public:
    util::Point d_x;
 
    util::Point d_xBegin;
 
    util::Point d_xa;
 
    double d_r;
 
    double d_l;
 
  public:
    Cylinder()
      : GeomObject("cylinder", ""),
        d_xBegin(util::Point()),
        d_xa(util::Point(1., 1., 1.)),
        d_r(0.),
        d_l(0.),
        d_x(util::Point()) {
    };
 
    Cylinder(double r, double l, util::Point x_begin, util::Point xa,
             std::string description = "")
      : GeomObject("cylinder", description),
        d_xBegin(x_begin),
        d_xa(xa / xa.length()),
        d_r(r),
        d_l(l),
        d_x(x_begin + 0.5 * l * (xa / xa.length())) {
    };
 
    Cylinder(double r, util::Point x_begin, util::Point xa,
              std::string description = "")
            : GeomObject("cylinder", description),
              d_xBegin(x_begin),
              d_xa(xa / xa.length()),
              d_r(r),
              d_l(xa.length()),
              d_x(x_begin + 0.5 * xa) {};
 
    Cylinder(const Cylinder &other)
      : GeomObject(other.d_name, other.d_description),
        d_xBegin(other.d_xBegin),
        d_xa(other.d_xa),
        d_r(other.d_r),
        d_l(other.d_l),
        d_x(other.d_x) {
      d_tags = other.d_tags;
    }
 
    Cylinder& operator=(const Cylinder& other) {
      if (this != &other) {
        // Copy base class members
        d_tags = other.d_tags;
        
        // Copy member variables
        d_xBegin = other.d_xBegin;
        d_xa = other.d_xa;
        d_r = other.d_r;
        d_l = other.d_l;
        d_x = other.d_x;
      }
      return *this;
    }
 
    void transform(const util::Point &translation, const double &scale, 
                  const double &angle, const util::Point &axis, const util::Point *rotationPoint) override {
      // Scale the dimensions
      d_r *= scale;
      d_l *= scale;
 
      const util::Point c0 = d_x;
      d_xBegin = mapSimilarity(d_xBegin, c0, translation, scale, angle, axis, rotationPoint);
      d_xa = util::rotate(d_xa, angle, axis);
      d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
    }
 
    double volume() const override;
 
    util::Point center() const override;
 
    std::pair<util::Point, util::Point> box() const override;
 
    std::pair<util::Point, util::Point>
    box(const double &tol) const override;
 
    double inscribedRadius() const override;
 
    double boundingRadius() const override;
 
    bool isInside(const util::Point &x) const override;
 
    bool isOutside(const util::Point &x) const override;
 
    bool isNear(const util::Point &x, const double &tol) const override;
 
    bool isNearBoundary(const util::Point &x, const double &tol,
                        const bool &within) const override;
 
    bool doesIntersect(const util::Point &x) const override;
 
    bool
    isInside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool
    isOutside(
      const std::pair<util::Point, util::Point> &box) const override;
 
    bool isNear(const std::pair<util::Point, util::Point> &box,
                const double &tol) const override;
 
    bool doesIntersect(
      const std::pair<util::Point, util::Point> &box) const override;
 
    std::string printStr(int nt, int lvl) const override;
 
    void print(int nt, int lvl) const override {
      std::cout << printStr(nt, lvl);
    };
 
    void print() const override { print(0, 0); };
  };
} // namespace geom
 
#endif // GEOM_GEOMOBJECTS_H