Homec4science
Diffusion CADD_mesher (master) Browse

rCADDMESH/src/common/point_container.hh
1c407f7780f6master

Last Change
Raw File

point_container.hh
View Options

/**
* @file point_container.hh
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Fri Apr 13 11:10:22 2012
*
* @brief
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __CADD_MESH_POINT_CONTAINER_HH__
#define __CADD_MESH_POINT_CONTAINER_HH__
#include <vector>
#include "string"
#include "sstream"
#include "fstream"
#include "common.hh"
#include "dumper_paraview.hh"
#include "dumper_lammps.hh"
#include <limits>
#include <map>
#include <array>
#include <stdexcept>
template <Uint DIM>
class Geometry;
/*! \class PointRef simple representation of a point.*/
template<Uint DIM>
class PointRef {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/*! Constructor
* \param _vals array of which the point is part
* \param _index position of the point in the array. E.g., for 3d points, the
n-th point would have index n-1, and \a not 3*(n-1)*/
PointRef(PointRef& other)
:is_light(other.is_light)
{
if (this->is_light) {
this->vals = other.vals;
this->index = other.index;
} else {
this->vals = new Real [DIM];
std::copy(other.vals, other.vals + DIM, this->vals);
this->index = 0;
}
};
PointRef(const PointRef& other)
:is_light(other.is_light)
{
this->vals = new Real [DIM];
std::copy(other.vals, other.vals + DIM, this->vals);
this->index = 0;
};
PointRef<DIM> & operator=(const PointRef<DIM>& other) {
if (this != &other) {
if (this->is_light) {
this->vals = new Real [DIM];
}
std::copy(other.vals, other.vals + DIM, this->vals);
this->is_light = false;
this->index = 0;
}
return *this;
}
PointRef<DIM> operator-(const PointRef<DIM>& other) {
PointRef<DIM> tmp;
std::copy(this->vals, this->vals + DIM, tmp.vals);
for (Uint i = 0 ; i < DIM ; ++i) {
tmp.vals[i] -= other[i];
}
return tmp;
}
PointRef<DIM> operator+(const PointRef<DIM>& other) {
PointRef<DIM> tmp;
std::copy(this->vals, this->vals + DIM, tmp.vals);
for (Uint i = 0 ; i < DIM ; ++i) {
tmp.vals[i] += other[i];
}
return tmp;
}
PointRef(Real * _vals=NULL, Uint _index=0)
: index(_index) {
if (_vals == NULL) {
this->is_light = false;
this->vals = new Real[DIM];
} else {
this->is_light = true;
this->vals = _vals + _index * DIM;
}
}
PointRef(std::initializer_list<Real> inlist)
:vals(new Real[DIM]), index(0), is_light(false)
{
if (inlist.size() != DIM) {
throw std::invalid_argument("Size of initializer list must be DIM");
}
std::copy(inlist.begin(), inlist.end(), this->vals);
}
PointRef(std::vector<Real> & coords)
:index(0), is_light(false) {
this->vals = new Real[DIM];
for (Uint i = 0; i < DIM ; ++i) {
this->vals[i] = coords.at(i);
}
}
PointRef copy() const {
PointRef out_pt;
out_pt = *this;
return out_pt;
}
bool operator==(const PointRef & other) const {
bool equal = true;
for (Uint i = 0 ; i < DIM ; ++i) {
equal &= (this->getComponent(i) == other.getComponent(i));
}
return equal;
}
virtual ~PointRef(){ if (!this->is_light) delete[] this->vals;};
Real norm() {return vector_norm<DIM>(this->vals);}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void make_me_independent() {
if (this-> is_light) {
Real new_vals[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
new_vals[i] = this->vals[i];
}
this->vals = new_vals;
this->is_light = false;
}
}
void print_coords(std::ostream & stream) const
{
std::ios_base::fmtflags current_flags = stream.flags();
stream << std::scientific << "(";
for (Uint i = 0; i < DIM -1; ++i) {
stream << this->vals[i] << ", ";
}
stream << this->vals[DIM -1] << ")";
stream.setf(current_flags);
}
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const
{
indent_space(stream, indent);
stream << "Point ";
stream.width(6);
stream << this->index << ": ";
this->print_coords(stream);
}
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/*! self-explanatory
*/
inline Real & operator[](Uint dir) {return this->vals[dir];}
inline const Real & operator[](Uint dir) const {return this->vals[dir];}
inline const Real & getComponent(Uint dir) const {return this->vals[dir];}
inline Real * getArray() {return this->vals;}
inline const Real * getArray() const {return this->vals;}
inline std::vector<Real> getVector() const {
return std::vector<Real>(this->vals, this->vals+DIM);
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
Real * vals;
Uint index;
bool is_light;
};
template<Uint DIM>
inline Real distance2(const PointRef<DIM> & pointA,
const PointRef<DIM> & pointB) {
Real dist2 = 0;
for (Uint i = 0; i < DIM ; ++i) {
Real manhattan = pointA.getComponent(i) - pointB.getComponent(i);
dist2 += manhattan * manhattan ;
}
return dist2;
}
template<Uint DIM>
inline Real distance(const PointRef<DIM> & pointA,
const PointRef<DIM> & pointB) {
return std::sqrt(distance2(pointA, pointB));
}
/*! \class PointContainer container for vectors
*/
template<Uint DIM>
class PointContainer {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor \param name by default ""
PointContainer(std::string _name = ""):name(_name), threshold(-1.),
modified(true){};
PointContainer(const PointContainer & other)
:name(other.name), values(other.values), threshold(-1.),
modified(true) {};
PointContainer(const std::vector<Real> & coords,
const std::string & _name = "")
:name(_name), values(coords), threshold(-1.), modified(true){};
virtual ~PointContainer(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
PointContainer & operator=(const PointContainer & other) {
this->name = other.name;
this->values = other.values;
this->modified = true;
return *this;
}
PointContainer operator+(const PointContainer & other) const {
PointContainer out_points;
for (Uint i = 0 ; i < this->getSize() ; ++i) {
out_points.addPoint(this->getConstPoint(i));
}
for (Uint i = 0 ; i < other.getSize() ; ++i) {
out_points.addPoint(other.getConstPoint(i));
}
return out_points;
}
PointContainer & operator+=(const PointContainer & other) {
for (Uint i = 0 ; i < other.getSize() ; ++i) {
this->addPoint(other.getConstPoint(i));
}
return *this;
}
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const
{
indent_space(stream, indent);
stream << "PointContainer '" << this->name << "' contains "
<< this->getSize() << " points:" << std::endl;
Uint size = this->getSize();
Uint max_size = 20;
Uint i;
for (i = 0; i < size && i < max_size; ++i) {
PointRef<DIM> printref = this->getConstPoint(i);
printref.printself(stream, indent + indent_incr);
stream << std::endl;
}
if (size > max_size) {
indent_space(stream, indent + indent_incr);
stream << "Point ..." << std::endl;
}
}
void computeBounds(Real * bounds) const {
for (Uint dir = 0 ; dir < DIM ; ++dir) {
bounds[2 * dir + 0] = std::numeric_limits<Real>::max();
bounds[2 * dir + 1] = -std::numeric_limits<Real>::max();
}
for (Uint point = 0 ; point < this->values.size() ; point+=DIM) {
for (Uint dir = 0 ; dir < DIM ; ++dir) {
keep_min(bounds[2*dir + 0], this->values.at(point+dir));
keep_max(bounds[2*dir + 1], this->values.at(point+dir));
}
}
}
std::vector<Real> computeBounds() const {
std::vector<Real> bounds;
bounds.resize(2*DIM);
this->computeBounds(&bounds[0]);
return bounds;
}
void filter(const Geometry<DIM> & geom, Real tol);
void dump_paraview(std::string filename="diagnostic_pointcontainer") const {
iohelper::DumperParaview dumper;
dumper.setMode(iohelper::TEXT);
dumper.setPoints(const_cast<Real*>(&this->values[0]),
DIM, this->values.size()/DIM, filename);
dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE);
dumper.setPrefix("./");
dumper.init();
try {
dumper.dump();
} catch (iohelper::IOHelperException & err) {
std::cout << err.what() << std::endl;
throw(err);
}
}
void dump_lammps(std::string filename, const Real & skin) {
Real bounds[2*DIM];
this->computeBounds(bounds);
for (Uint i = 0 ; i < DIM ; ++i) {
bounds[2*i ] -= skin;
bounds[2*i + 1] += skin;
}
this->dump_lammps(filename, bounds);
}
void dump_lammps(std::string filename, std::vector<Real> bounds) const{
this->dump_lammps(filename, &bounds[0]);
}
void dump_lammps(std::string filename="diagnostic_pointcontainer",
Real * set_bounds = NULL) const {
Real bounds[2*DIM];
if (set_bounds == NULL) {
this->computeBounds(bounds);
} else {
for (Uint i = 0 ; i < DIM*2 ; ++i) {
bounds[i] = set_bounds[i];
}
}
iohelper::DumperLammps<iohelper::atomic> dumper(bounds);
dumper.setPoints(const_cast<Real*>(&this->values[0]),
DIM, this->values.size()/DIM, filename);
dumper.setPrefix("./");
dumper.setConnectivity(NULL, iohelper::POINT_SET, 0, iohelper::C_MODE);
dumper.init();
dumper.dump();
}
void dump_text(std::string filename="diagnostic_pointcontainer",
std::string commentary="") const {
std::ofstream ofile(filename.c_str());
if (ofile.good()) {
ofile << "# file: " << filename <<std::endl;
if (commentary.length()>0) {
ofile << "# " << commentary << std::endl;
}
}
ofile << std::setprecision(16) << std::scientific << std::showpos;
Uint current_point = 0;
for (Uint i = 0 ; i < this->getSize(); ++i) {
for (Uint j = 0 ; j < DIM - 1 ; ++j, ++current_point) {
ofile << this->values.at(current_point) << "\t";
}
ofile << this->values.at(current_point++);
ofile << std::endl;
}
ofile.close();
}
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/*! \param coords array of coordinates to add
*/
void addPoint(const Real * coords)
{
for (Uint i = 0; i < DIM; ++i) {
this->values.push_back(coords[i]);
}
this->modified = true;
}
void addPoint(std::vector<Real> & coords) throw (std::string) {
if (coords.size() < DIM) {
std::stringstream error_msg;
error_msg << "Size of coords should be " << DIM
<< " but it's " << coords.size();
throw error_msg.str();
}
this->addPoint(&coords[0]);
}
inline void addPoint(const PointRef<DIM> & point) {
this->addPoint(&point.getComponent(0));}
/// returns the number of points in the container
inline Uint getSize() const {return this->values.size()/DIM;}
/// returns the container name
std::string getName() const {return this->name;}
/// sets the container name
void setName(std::string _name) {this->name = _name;}
/*! returns a PointRef to the point at position \a index. Throws a string if
* the index is out of bounds
* \param index */
PointRef<DIM> getPoint(Uint index) throw (std::runtime_error)
{
if (index >= this->getSize()) {
std::stringstream error_msg;
error_msg << "index " << index << " is out of bounds. Size of vector "
<< this->name << " is only " << this->getSize();
throw std::runtime_error(error_msg.str());
}
return PointRef<DIM>(&this->values[0], index);
}
const PointRef<DIM> getConstPoint(Uint index) const throw (std::runtime_error)
{
if (index >= this->getSize()) {
std::stringstream error_msg;
error_msg << "index " << index << " is out of bounds. Size of vector "
<< this->name << " is only " << this->getSize();
throw std::runtime_error(error_msg.str());
}
const PointRef<DIM> ret_ref(const_cast<Real *>(&values[0]), index);
return ret_ref;
}
const PointRef<DIM> operator[](Uint index) const throw (std::string) {
return this->getConstPoint(index);
}
PointRef<DIM> operator[](Uint index) throw (std::string) {
return this->getConstPoint(index);
}
/*! returnd the vector of raw data. You'll probably never use this.*/
const std::vector<Real> & getVector() const {return this-> values;}
const Real * getArray() const {return &this->values[0];}
/*! only use this if you know what you're doing. This can fuck up values
and thus render the container useless */
std::vector<Real> & getNonConstVector() {return this->values;}
void deletePoint(Uint index) {
Uint nb_points = this->getSize();
for (Uint i = 0 ; i < DIM ; ++i) {
this->values[DIM*index+i] = this->values[DIM*(nb_points-1) + i];
}
for (Uint i = 0 ; i < DIM ; ++i) {
this->values.pop_back();
}
this->modified = true;
}
void clear();
Uint reportDuplicates(Real threshold);
Uint cleanDuplicates(Real threshold);
// check whether exactly this point (no rounding error) is in the container
bool containsExacly(const PointRef<DIM> &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
template<bool do_delete> Uint sniffDuplicates();
void createSearchGrid();
void setThreshold(Real thresh);
std::string name;
std::vector<Real> values;
// search grid related members
Real threshold;
std::map<std::array<int, DIM>, std::vector<Uint> > grid;
std::map<Uint, std::array<int, DIM> > reverse_grid;
std::vector<std::array<int, DIM> > boxes;
bool modified; //grid needs to be recomputed if modifications happened
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
/// standard output stream operator
template<Uint DIM>
inline std::ostream & operator <<(std::ostream & stream, const PointContainer<DIM> & _this)
{
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
/// standard output stream operator
template<Uint DIM>
inline std::ostream & operator <<(std::ostream & stream, const PointRef<DIM> & _this)
{
_this.printself(stream);
return stream;
}
#endif /* __CADD_MESH_POINT_CONTAINER_HH__ */
c4science · Help