geometry.h

Go to the documentation of this file.

#include "palabos3D.h"
#include "palabos3D.hh"
#include <cmath>
#include <limits>
 
#ifndef GEOMETRY_H
#define GEOMETRY_H
 
namespace geom {
 
class BoxDomain : public plb::DomainFunctional3D {
public:
  explicit BoxDomain(plb::Box3D box) : box{box} {}
 
  bool operator()(plb::plint x, plb::plint y, plb::plint z) const override {
    return plb::contained(x, y, z, box);
  }
 
  const plb::Box3D& bounding_box() const { return box; }
 
  BoxDomain *clone() const override { return new BoxDomain(*this); }
private:
  plb::Box3D box;
};
 
class EllipseDomain : public plb::DomainFunctional3D {
public:
  EllipseDomain(double x, double y, double rx, double ry)
      : EllipseDomain(x, y, 0, rx, ry, std::numeric_limits<double>::max())
  {}
 
  EllipseDomain(double cx, double cy, double cz, double rx, double ry, double rz)
      : cx{ cx },
        cy{ cy },
        cz{ cz },
        rx{ rx },
        ry{ ry },
        rz{ rz }
  {}
 
  bool operator()(plb::plint x, plb::plint y, plb::plint z) const override {
    auto dx = x - cx;
    auto dy = y - cy;
    auto dz = z - cz;
    return (dx*dx)/(rx*rx) + (dy*dy)/(ry*ry) + (dz*dz)/(rz*rz) < 1;
  }
 
  EllipseDomain *clone() const override { return new EllipseDomain(*this); }
 
private:
  double cx = 0;
  double cy = 0;
  double cz = 0;
  double rx = 0;
  double ry = 0;
  double rz = 0;
};
 
class Union : public plb::DomainFunctional3D {
public:
  Union(std::initializer_list<std::reference_wrapper<plb::DomainFunctional3D>>
            domains)
      : domains{domains}
  {}
 
  bool operator()(plb::plint x, plb::plint y, plb::plint z) const override {
    auto predicate = [x, y, z](const plb::DomainFunctional3D &domain) {
      return domain(x, y, z);
    };
    return std::any_of(domains.begin(), domains.end(), predicate);
  }
 
  Union *clone() const override { return new Union(*this); };
 
private:
  std::vector<std::reference_wrapper<plb::DomainFunctional3D>> domains;
};
 
class Difference : public plb::DomainFunctional3D {
public:
  Difference(
      std::initializer_list<std::reference_wrapper<plb::DomainFunctional3D>> domains)
      : domains{ domains } {}
 
  bool operator()(plb::plint x, plb::plint y, plb::plint z) const override {
    if (domains.empty())
      return false;
 
    if (!domains.front()(x, y, z))
      // No need to subtract subsequent domains if the first domain is empty.
      return false;
 
    auto predicate = [x, y, z](const plb::DomainFunctional3D& domain) {
      return domain(x, y, z);
    };
    return std::none_of(std::next(domains.begin()), domains.end(), predicate);
  }
 
  Difference *clone() const override { return new Difference(*this); }
 
private:
  std::vector<std::reference_wrapper<plb::DomainFunctional3D>> domains;
};
 
class Intersection : public plb::DomainFunctional3D {
public:
  Intersection(
      std::initializer_list<std::reference_wrapper<plb::DomainFunctional3D>>
          domains)
      : domains{domains} {}
 
  bool operator()(plb::plint x, plb::plint y, plb::plint z) const override {
    if (domains.empty())
      // No domains have an empty intersection, whereas std::all_of returns
      // true for emtpy ranges.
      return false;
 
    auto predicate = [x, y, z](const plb::DomainFunctional3D &domain) {
      return domain(x, y, z);
    };
    return std::all_of(domains.begin(), domains.end(), predicate);
  }
 
  Intersection *clone() const override { return new Intersection(*this); }
 
private:
  std::vector<std::reference_wrapper<plb::DomainFunctional3D>> domains;
};
 
} // namespace geom
 
#endif /* GEOMETRY_H */