Page MenuHomec4science

point_container.hh
No OneTemporary

File Metadata

Created
Sat, Aug 17, 18:53

point_container.hh

/**
* @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>
/*! \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(const 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<DIM> & operator=(const PointRef<DIM>& other) {
if (this != &other) {
this->vals = new Real [DIM];
std::copy(other.vals, other.vals + DIM, this->vals);
this->is_light = false;
this->index = 0;
}
return *this;
}
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::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);
}
}
bool operator==(const PointRef & other) {
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;};
/* ------------------------------------------------------------------------ */
/* 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;
}
}
/// 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 << ": ";
stream << std::scientific << "(";
for (Uint i = 0; i < DIM -1; ++i) {
stream << this->vals[i] << ", ";
}
stream << this->vals[DIM -1] << ")" << std::endl;
}
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/*! self-explanatory
*/
inline const Real & getComponent(Uint index) const {return this->vals[index];}
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>
class Geometry;
/*! \class PointContainer container for vectors
*/
template<Uint DIM>
class PointContainer {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor \param name by default ""
PointContainer(std::string _name = ""):name(_name){};
PointContainer(const PointContainer & other)
:name(other.name), values(other.values) {};
virtual ~PointContainer(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
PointContainer & operator=(const PointContainer & other) {
this->name = other.name;
this->values = other.values;
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);
}
if (size > max_size) {
indent_space(stream, indent + indent_incr);
stream << "Point ..." << std::endl;
}
}
void computeBounds(Real * bounds) {
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() {
std::vector<Real> bounds;
bounds.resize(2*DIM);
this->computeBounds(&bounds[0]);
return bounds;
}
void filter(const Geometry<DIM> & geom) {
PointContainer temp;
for (Uint point = 0 ; point < this->getSize() ; ++point) {
PointRef<DIM> current_point = this->getPoint(point);
if (geom.is_inside(current_point)) {
temp.addPoint(current_point);
}
}
this->values = std::move(temp.values);
}
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, std::vector<Real> bounds) {
this->dump_lammps(filename, &bounds[0]);
}
void dump_lammps(std::string filename="diagnostic_pointcontainer",
Real * set_bounds = NULL) {
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]);
}
}
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::string)
{
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 error_msg.str();
}
return PointRef<DIM>(&this->values[0], index);
}
const PointRef<DIM> getConstPoint(Uint index) const throw (std::string)
{
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 error_msg.str();
}
const PointRef<DIM> ret_ref(const_cast<Real *>(&values[0]), index);
return ret_ref;
}
/*! returnd the vector of raw data. You'll probably never use this.*/
const std::vector<Real> & getVector() const {return this-> values;}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
std::string name;
std::vector<Real> values;
};
/* -------------------------------------------------------------------------- */
/* 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__ */

Event Timeline