ReadQuake.h

Go to the documentation of this file.

#ifndef ReadQuake_HEADER
#define ReadQuake_HEADER

#include <assert.h>
#include <vector>
#include <geometry/Point.h>

class QuakeReader {
  public:
    QuakeReader(const char *filename);
    ~QuakeReader();

    typedef size_t Dimension;
    typedef size_t ID;

    struct point {
      ID id;
      std::vector<double> coords;
      std::vector<double> attributes;
      bool isBdry;
    };

    struct element {
      ID id;
      Dimension dim;
      std::vector<ID> subElements; // IDs of dim-1 elements
      std::vector<std::pair<Dimension,ID> > orphans; // <dim, id> of lower-d elements
      bool isBdry;
    };

    struct point_iterator {
      void operator ++ () {
        id_++;
        assert(id_ <= source_->afterLastID());
      }
      int operator- (const point_iterator& other) const {
        return id_ - other.id_;
      }
      bool operator!= (const point_iterator& other) const {
        return id_ != other.id_;
      }
      point operator*() const {
        return source_->get_point(id_);
      }

      protected:
      friend class QuakeReader;
      point_iterator(const QuakeReader *r, ID ii)
        : source_(r), id_(ii) { }
      const QuakeReader *source_;
      ID id_;
    };

    typedef std::vector<element>::const_iterator element_iterator;

    Dimension dim() const { return dim_; }
    bool hasConstraints() const;

    size_t npoints() const { return npoints_; }
    ID firstID() const { return firstID_; }
    ID afterLastID() const { return firstID_ + npoints(); }
    point_iterator begin_points() const;
    point_iterator end_points() const;
    point_iterator try_get_point(ID) const;
    point get_point(ID) const;
    // Dimension dim() const
    const double *get_x(ID) const;
    size_t nPointAttrs() const { return nattrs_; }
    const double *get_attrs(ID) const;
    bool isBdry(ID) const;

    template <size_t d> struct PointIterator : point_iterator {
      typedef point_iterator super;
      PointIterator(const point_iterator& it): super(it) { }
      void operator ++ () { super::operator++(); }
      int operator- (const PointIterator<d>& other) const {
        return super::operator-(other);
      }
      bool operator != (const PointIterator<d>& other) const {
        return super::operator!=(other);
      }
      Geometry::Point<d> operator*() const {
        return Geometry::Point<d>(source_->get_x(id_), d);
      }
    };

    template <size_t d> PointIterator<d> beginPoints() const {
      return PointIterator<d>(begin_points());
    }
    template <size_t d> PointIterator<d> endPoints() const {
      return PointIterator<d>(end_points());
    }

    size_t nelements(Dimension) const;
    element_iterator begin_elements(Dimension) const;
    element_iterator end_elements(Dimension) const;
    element_iterator try_get_element(std::pair<Dimension,ID>) const;
    const element& get_element(std::pair<Dimension,ID> p) const;

  private:
    Dimension dim_;

    /* We're very careful with the allocation here, because it turns out the
     * algorithm is so fast that I/O matters to the runtime. */
    ID firstID_;
    size_t npoints_;
    size_t nattrs_;
    double *pointX_;
    double *pointAttrs_;
    bool *pointIsBdry_;

    std::vector<std::vector<element> > elements_;

    struct NodeHeader {
      int n;
      int d;
      int nattr;
      int bdry;
    };
    void allocatePoints(const NodeHeader&);
};

#endif

---