diff --git a/src/common/common.hh b/src/common/common.hh
index 1f60f28..6492663 100644
--- a/src/common/common.hh
+++ b/src/common/common.hh
@@ -1,192 +1,192 @@
/**
* @file common.hh
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Wed Apr 4 15:59:37 2012
*
* @brief common parts for the cadd_mesher
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __cadd_mesh_COMMON_HH__
#define __cadd_mesh_COMMON_HH__
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
typedef double Real;
typedef unsigned int Uint;
const static Uint twoD = 2;
const static Uint threeD = 3;
const static std::string default_flags = std::string("qhull d Qt");
const static Uint indent_incr = 2;
void indent_space(std::ostream & stream, int indent = 0);
template <typename numeric>
inline void keep_max(numeric & value, const numeric & comparison) {
if (comparison > value) {
value = comparison;
}
}
template <typename numeric>
inline void keep_min(numeric & value, const numeric & comparison) {
if (comparison < value) {
value = comparison;
}
}
template <typename numeric>
inline numeric max(const numeric & x, const numeric & y) {
if (x > y) return x;
return y;
}
template <typename numeric>
inline numeric min(const numeric & x, const numeric & y) {
if (x > y) return y;
return x;
}
template <typename numeric>
inline numeric square(numeric val) {return val*val;}
template <typename el_typo>
void print_vector(const std::vector<el_typo> & vec, const std::string & name,
Uint stride = 1, int indent = 0, std::ostream & stream = std::cout){
indent_space(stream, indent);
stream << "Vector '" << name << "': " <<std::endl;
const Uint max_size = 12;
Uint size = vec.size()/stride;
Uint i;
for (i = 0; i < size && i < max_size; ++i) {
stream.width(7);
indent_space(stream, indent+indent_incr);
stream << i << " : " ;
for (Uint j = 0; j < stride-1 ; ++j) {
stream.width(4);
stream << std::showpos << vec.at(stride*i +j) << ",\t" ;
}
stream.width(4);
stream << vec.at(stride*i+stride-1) << std::endl;
}
if (size > max_size) {
stream.width(7);
indent_space(stream, indent+indent_incr);
stream << i << " : ";
for (Uint j = 0; j < stride-1; ++j) {
stream.width(4);
stream << std::showpos << "..." << ",\t";
}
stream.width(4);
stream << "..." << std::endl;
}
}
template <typename el_typo>
void print_array(el_typo * array, Uint length, const std::string & name, std::ostream & stream=std::cout) {
stream << "Array '" << name << "': (" ;
for (Uint i = 0; i < length-1 ; ++i) {
stream << array[i] << ", ";
}
stream << array[length-1] << ")";
}
//This function sets a vec1's components to vec2's components
template <Uint DIM>
inline void vector_set(Real * vec1, const Real * vec2) {
for (Uint i = 0; i < DIM ; ++i) {
vec1[i] = vec2[i];
}
}
//This function computes the norm of a vector
template <Uint DIM>
inline Real vector_norm(const Real vec[DIM]) {
Real norm2 = 0;
for (Uint i = 0; i < DIM ; ++i) {
norm2 += square(vec[i]);
}
return std::sqrt(norm2);
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
inline Real dot_prod(const Real vec1[DIM], const Real vec2[DIM]) {
Real prod = 0;
for (Uint i = 0 ; i < DIM ; ++i) {
prod += vec1[i]*vec2[i];
}
return prod;
}
template<Uint DIM>
inline Real norm_cross_prod(const Real vec1[DIM], const Real vec2[DIM]) {
Real prod = 0;
prod += vec1[0]*vec2[1] - vec1[1]*vec2[0];
if (DIM == twoD) {
return std::abs(prod);
}
- if (DIM == 3) {
+ if (DIM == threeD) {
prod *= prod;
prod += square(vec1[1]*vec2[2] - vec1[2]*vec2[1]);
prod += square(vec1[2]*vec2[0] - vec1[0]*vec2[2]);
}
return std::sqrt(prod);
}
// stretches a vector
template<Uint DIM>
inline void vector_stretch(Real * vector, const Real * factors) {
for (Uint i = 0; i < DIM ; ++i) {
vector[i] *= factors[i];
}
}
template<Uint DIM>
inline void vector_scale(Real * vector, const Real & factor) {
for (Uint i = 0; i < DIM ; ++i) {
vector[i] *= factor;
}
}
//This function adds factor times addition to vec1
template <Uint DIM>
inline void vector_incr(Real * vec1, const Real * addition, const Real & factor=1) {
for (Uint i = 0; i < DIM ; ++i) {
vec1[i] += factor * addition[i];
}
}
// vector from point 1 to point 2
template<Uint DIM>
inline void vector_difference(Real * vec, const Real * from, const Real * to,
const Real & factor = 1) {
vector_set<DIM>(vec, to);
vector_incr<DIM>(vec, from, -1);
vector_scale<DIM>(vec, factor);
}
// cross product
inline void cross_prod(const Real in0[threeD], const Real in1[threeD],
Real out[threeD]) {
out[0] = in0[1]*in1[2] - in0[2]*in1[1];
out[1] = in0[2]*in1[0] - in0[0]*in1[2];
out[2] = in0[0]*in1[1] - in0[1]*in1[0];
}
#endif /* __cadd_mesh_COMMON_HH__ */
diff --git a/src/common/point_container.hh b/src/common/point_container.hh
index 5bdf445..b7bc0ee 100644
--- a/src/common/point_container.hh
+++ b/src/common/point_container.hh
@@ -1,448 +1,448 @@
/**
* @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>
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(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) {
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> 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::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) 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);
stream << std::endl;
}
/* ------------------------------------------------------------------------ */
/* 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;
}
/*! \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;
}
/// 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) 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, 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::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;
}
const PointRef<DIM> operator[](Uint index) const 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;}
/*! 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__ */
diff --git a/src/geometry/block.cc b/src/geometry/block.cc
index f1cbf9d..e9a3df8 100644
--- a/src/geometry/block.cc
+++ b/src/geometry/block.cc
@@ -1,367 +1,367 @@
#include "block.hh"
#include "block_surface_manipulations.hh"
#include <limits>
#include <cmath>
#include <vector>
#include <deque>
#include <utility>
#include <tuple>
#include <algorithm>
#include <cadd_mesh.hh>
#include <exception>
#include <sstream>
class BlockException: public std::exception {
public:
BlockException(std::stringstream & _what):what_stream(_what.str()) {};
BlockException(std::string &_what):what_stream(_what) {};
~BlockException() throw() {};
virtual const char* what() const throw() {
return this->what_stream.c_str();
}
private:
std::string what_stream;
};
template<Uint DIM>
Block<DIM>::Block(const std::vector<Real> & _min_max, std::string _name)
:Geometry<DIM>(_name) {
if (_min_max.size() != 2*DIM) {
std::stringstream error;
error << "The length of 'min_max' should be = " << 2*DIM;
throw BlockException(error);
}
for (Uint i = 0 ; i < 2*DIM ; ++i) {
this->min_max[i] = _min_max[i];
}
}
template<Uint DIM>
Block<DIM>::Block(const Real * min_max, std::string _name)
:Geometry<DIM>(_name) {
for (Uint i = 0; i < 2 * DIM; ++i) {
this->min_max[i] = min_max[i];
}
}
template<Uint DIM>
Block<DIM>::Block(const Block & other)
:Geometry<DIM>(other.name) {
for (Uint i = 0 ; i < 2 * DIM ; ++i) {
this->min_max[i] = other.min_max[i];
}
}
template <Uint DIM>
void Block<DIM>::shift(const PointRef<DIM> & offset) {
for (Uint i = 0 ; i < DIM ; ++i) {
this->min_max[2 * i + 0] += offset.getComponent(i);
this->min_max[2 * i + 1] += offset.getComponent(i);
}
}
template<Uint DIM>
InsideObject Block<DIM>::from_border(const Real * point, Real tol) const {
Real lower_violation = -std::numeric_limits<Real>::max();
Real upper_violation = -std::numeric_limits<Real>::max();
for (Uint i = 0; i < DIM; ++i) {
keep_max(lower_violation, this->min_max[2*i] - point[i]);
keep_max(upper_violation, point[i] - this->min_max[2*i+1]);
}
return {max(lower_violation, upper_violation),
max(lower_violation-tol, upper_violation+tol) <= 0};
}
template<Uint DIM, bool eval_min, bool eval_max>
inline void in_direction_helper(Real & min_val, Real & max_val,
const Real * direction, Real * point) {
Real dot_prod = 0;
Real norm = 0;
for (Uint i = 0; i < DIM; ++i) {
dot_prod += point[i] * direction[i];
norm += direction[i] * direction[i];
}
dot_prod /= sqrt(norm);
if (eval_min) {
keep_min(min_val, dot_prod);
}
if (eval_max) {
keep_max(max_val, dot_prod);
}
}
template<Uint DIM>
Real Block<DIM>::min_in_direction(const Real * dir, const Real & unit) const {
Real point[DIM];
Real min_val = std::numeric_limits<Real>::max();
for (Uint i = 0; i < 2; ++i) {
point[0] = min_max[i];
for (Uint j = 0; j < 2; ++j) {
point[1] = min_max[2 + j];
if (DIM == 3) {
for (Uint k = 0; k < 2; ++k) {
point[2] = min_max[4 + k];
in_direction_helper<DIM, true, false>(min_val, min_val, dir, point);
}
} else if (DIM == 2) {
in_direction_helper<DIM, true, false>(min_val, min_val, dir, point);
}
}
}
return min_val/unit;
}
template<Uint DIM>
Real Block<DIM>::max_in_direction(const Real * dir, const Real & unit) const {
Real point[DIM];
Real max_val = -std::numeric_limits<Real>::max();
for (Uint i = 0; i < 2; ++i) {
point[0] = min_max[i];
for (Uint j = 0; j < 2; ++j) {
point[1] = min_max[2 + j];
if (DIM == 3) {
for (Uint k = 0; k < 2; ++k) {
point[2] = min_max[4 + k];
in_direction_helper<DIM, false, true>(max_val, max_val, dir, point);
}
} else if (DIM == 2) {
in_direction_helper<DIM, false, true>(max_val, max_val, dir, point);
}
}
}
return max_val/unit;
}
template<Uint DIM>
Real Block<DIM>::size_in_direction(const Real * dir, const Real & unit) const {
Real point[DIM];
Real min_val = std::numeric_limits<Real>::max();
Real max_val = -std::numeric_limits<Real>::max();
for (Uint i = 0; i < 2; ++i) {
point[0] = min_max[i];
for (Uint j = 0; j < 2; ++j) {
point[1] = min_max[2 + j];
if (DIM == 3) {
for (Uint k = 0; k < 2; ++k) {
point[2] = min_max[4 + k];
in_direction_helper<DIM, true, true>(min_val, max_val, dir, point);
}
} else if (DIM == 2) {
in_direction_helper<DIM, true, true>(min_val, max_val, dir, point);
}
}
}
return (max_val - min_val)/unit;
}
template<Uint DIM>
void Block<DIM>::printself(std::ostream & stream, int indent) const {
indent_space(stream, indent);
stream << "Block " << this->get_name() << ": " << std::endl;
indent_space(stream, indent);
stream << "x_min " << this->min_max[0] << ", x_max = " << this->min_max[1] << std::endl;
indent_space(stream, indent);
stream << "y_min " << this->min_max[2] << ", y_max = " << this->min_max[3] << std::endl;
if (DIM == 3) {
indent_space(stream, indent);
stream << "z_min " << this->min_max[4] << ", z_max = " << this->min_max[5] << std::endl;
}
}
template<Uint DIM>
inline void enforce_slave(const PointRef<DIM> & point,
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real & step_size,
PointContainer<DIM> & points) {
Real ref_sq = step_size * auxiliary.interpolate(point, refinement);
ref_sq *= ref_sq;
std::deque<Uint> delete_points;
for (Uint point_id = 0 ; point_id < points.getSize() ; ++point_id) {
if (ref_sq > distance2(point, points.getPoint(point_id))) {
delete_points.push_front(point_id);
}
}
auto point_it = delete_points.begin();
const auto point_end = delete_points.end();
for(; point_it != point_end; ++point_it) {
points.deletePoint(*point_it);
}
points.addPoint(point);
}
template <Uint DIM>
inline bool masterdom(const PointRef<DIM> & point, const Real * min_max,
const bool periodicity[DIM], bool master[DIM], bool & slave) {
bool is_master = false;
bool is_not_master = false;
slave = false;
for (Uint i = 0 ; i < DIM ; ++i) {
master[i] = (point[i] == min_max[2*i]) && periodicity[i];
// not a master if i'm at the max of any periodic direction
is_not_master |= (point[i] == min_max[2*i + 1]) && periodicity[i];
is_master |= master[i];
slave |= (point[i] == min_max[2*i + 1]) && periodicity[i];
}
return is_master && !is_not_master;
}
template <Uint DIM>
inline void drive_slaves(const PointRef<DIM> & point, const Real * min_max,
const bool periodicity[DIM],
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real & step_size,
PointContainer<DIM> & points) {
bool master[DIM];
bool slave;
if (masterdom(point, min_max, periodicity, master, slave)) {
PointContainer<DIM> slaves;
slaves.addPoint(point);
PointRef<DIM> tmp;
for (Uint i = 0 ; i < DIM ; ++i) {
if ( master[i]) {
Uint current_nb_slaves = slaves.getSize();
for (Uint slave_num = 0 ; slave_num < current_nb_slaves ; ++slave_num) {
tmp = slaves.getConstPoint(slave_num);
tmp[i] += min_max[2*i+1]-min_max[2*i];
slaves.addPoint(tmp);
}
}
}
for (Uint slave_num = 0 ; slave_num < slaves.getSize() ; ++slave_num) {
enforce_slave(slaves.getConstPoint(slave_num),
auxiliary,
refinement,
step_size,
points);
}
} else if (!slave) {
points.addPoint(point);
}
}
template <Uint DIM>
inline void micro_add_surface_point(PointRef<DIM> & point_hi,
PointRef<DIM> & point_lo,
const Uint xi,
const Real * min_max,
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const bool periodicity[DIM],
const Real & repr_const,
PointContainer<DIM> & points) {
point_hi[xi] = point_lo[xi] = min_max[2*xi];
if (check_if_ok(point_lo, auxiliary, refinement, repr_const, points)) {
drive_slaves(point_lo, min_max, periodicity, auxiliary, refinement,
repr_const, points);
}
if (check_if_ok(point_hi, auxiliary, refinement, repr_const, points)) {
drive_slaves(point_hi, min_max, periodicity, auxiliary, refinement,
repr_const, points);
}
point_hi[xi] = point_lo[xi] = min_max[2*xi+1];
if (!periodicity[xi]) {
add_if_ok(point_hi, auxiliary, refinement, repr_const, points);
add_if_ok(point_lo, auxiliary, refinement, repr_const, points);
}
}
template <Uint DIM>
void recurse_through_dimensions(PointRef<DIM> & point,
const Uint xi,
const Uint my_offset_dim,
const Real * min_max,
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const bool periodicity[DIM],
const Real step_size[DIM],
const Real & repr_const,
PointContainer<DIM> & points,
Uint * starts) {
Uint my_dim = (xi+my_offset_dim) % DIM;
Real length = min_max[2*my_dim + 1] - min_max[2*my_dim];
Uint nb_steps = length/step_size[my_dim];
PointRef<DIM> point_lo = point;
PointRef<DIM> point_hi = point;
// We mesh from the outside inwards, to make sure that corners/edges get
// points
for (Uint step = starts[my_offset_dim-1]; step < nb_steps/2+1; ++step) {
point_lo[my_dim] = min_max[2*my_dim] + step * step_size[my_dim];
point_hi[my_dim] = min_max[2*my_dim + 1] - step * step_size[my_dim];
if (DIM == my_offset_dim + 1) {
micro_add_surface_point(point_hi, point_lo, xi, min_max, auxiliary,
refinement, periodicity, repr_const, points);
} else {
recurse_through_dimensions(point_lo, xi, my_offset_dim + 1, min_max,
auxiliary, refinement, periodicity, step_size,
repr_const, points, starts);
recurse_through_dimensions(point_hi, xi, my_offset_dim + 1, min_max,
auxiliary, refinement, periodicity, step_size,
repr_const, points, starts);
}
}
}
template <Uint DIM>
void Block<DIM>::generate_surface_points(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real step_size[DIM],
const Real & repr_const,
const bool periodicity[DIM],
PointContainer<DIM> & points) const {
PointRef<DIM> point;
Uint starts[DIM] = {0};
for (Uint xi = 0 ; xi < DIM ; ++xi) {
starts[DIM-1-xi] = 1;
recurse_through_dimensions(point, xi, 1, this->min_max, auxiliary,
refinement, periodicity, step_size,
repr_const, points, starts);
}
points.cleanDuplicates(1.);
}
template<Uint DIM>
void Block<DIM>::complement_periodically(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
PointContainer<DIM> & points,
int direction,
Uint dim) const {
// relevant coordinate of the good and the bad side
Real ok_coord, bad_coord;
if (direction > 0) {
ok_coord = this->min_max[2*dim + 0];
bad_coord = this->min_max[2*dim + 1];
} else {
ok_coord = this->min_max[2*dim + 1];
bad_coord = this->min_max[2*dim + 0];
}
Real tol = std::max(fabs(ok_coord), 1.);
tol = 2000*tol*std::numeric_limits<Real>::epsilon();
PointRef<DIM> point;
for (Uint i = 0; i < points.getSize() ; ++i) {
point = points.getPoint(i);
if (fabs(point.getComponent(dim) - ok_coord) < tol) {
point.getArray()[dim] = bad_coord;
Real local_refinement = auxiliary.interpolate(point, refinement);
if (local_refinement <= 1.) {
points.addPoint(point);
}
}
}
}
-template class Block<2>;
-template class Block<3>;
+template class Block<twoD>;
+template class Block<threeD>;
diff --git a/src/geometry/cylinder.cc b/src/geometry/cylinder.cc
index 046c870..0bcd01a 100644
--- a/src/geometry/cylinder.cc
+++ b/src/geometry/cylinder.cc
@@ -1,233 +1,262 @@
#include "cylinder.hh"
#include <exception>
#include <sstream>
#include <limits>
/* -------------------------------------------------------------------------- */
class CylinderException: public std::exception {
public:
CylinderException(std::stringstream & _what): what_stream(_what.str()) {};
CylinderException(std::string & _what): what_stream(_what) {};
~CylinderException() throw() {};
virtual const char* what() const throw() {
return this->what_stream.c_str();
}
private:
std::string what_stream;
};
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Cylinder<DIM>::Cylinder(const std::vector<Real> & centre, const Real & radius_,
const std::vector<Real> & axis_, std::string name)
:Geometry<DIM>(name)
{
std::stringstream error;
if (DIM < 3) {
error << "This constructor only makes sense for more than 2 dimensions.";
throw CylinderException(error);
}
Uint size_centre = centre.size();
if (size_centre != DIM) {
error << "The length of centre should have been " << DIM
<< " but I got " << size_centre << "." ;
throw CylinderException(error);
}
Uint size_axis = axis_.size();
if (size_axis != DIM) {
error << "The length of axis should have been " << DIM
<< " but I got " << size_centre << "." ;
throw CylinderException(error);
}
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->base[dir] = centre[dir];
this->axis[dir] = axis_[dir];
}
this->radius = radius_;
this->length = vector_norm<DIM>(this->axis);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Cylinder<DIM>::Cylinder(const std::vector<Real> & centre, const Real & radius_,
std::string name)
:Geometry<DIM>(name)
{
std::stringstream error;
if (DIM > 2) {
error << "This constructor only makes sense for less than 3 dimesions.";
throw CylinderException(error);
}
Uint size_centre = centre.size();
if (size_centre != DIM) {
error << "The length of centre should have been " << DIM
<< " but I got " << size_centre << "." ;
throw CylinderException(error);
}
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->base[dir] = centre[dir];
this->axis[dir] = 0;
}
this->radius = radius_;
this->length = vector_norm<DIM>(this->axis);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Cylinder<DIM>::Cylinder(const PointRef<DIM> & centre, const Real & radius,
const Real * axis_, std::string name)
:Geometry<DIM>(name)
{
this->base = centre;
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->axis[dir] = axis_[dir];
}
this->radius = radius;
this->length = vector_norm<DIM>(this->axis);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Cylinder<DIM>::~Cylinder(){};
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Cylinder<DIM>::printself(std::ostream & stream, int indent) const {
indent_space(stream, indent);
stream << "Cylinder " << this->get_name() << ": " << std::endl;
indent_space(stream, indent);
stream << "base : " << this->base << std::endl;
indent_space(stream, indent);
- stream << "axis : " << PointRef<DIM>(this->axis) << std::endl;
+ stream << "axis : " << PointRef<DIM>(const_cast<Real*>(this->axis))
+ << std::endl;
indent_space(stream, indent);
stream << "radius : " << this->radius;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Cylinder<DIM>::shift(const PointRef<DIM> & offset) {
for (Uint i = 0 ; i < DIM ; ++i) {
this->base[i] += offset[i];
}
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
InsideObject Cylinder<DIM>::from_border(const Real * point, Real tol) const {
- PointRef<DIM> relpos = PointRef<DIM>(point) - this->base;
- Real violation = norm_cross_prod<DIM>(relpos, axis) - this->radius;
+ PointRef<DIM> relpos = PointRef<DIM>(const_cast<Real *>(point)) - this->base;
+ Real violation = norm_cross_prod<DIM>(relpos.getArray(), axis) - this->radius;
if (DIM == threeD) {
- Real dot_p = dot_prod(relpos, axis);
+ Real dot_p = dot_prod<DIM>(relpos.getArray(), axis);
keep_max(violation, -dot_p);
- keep_max(violation, dot_p - this-length);
+ keep_max(violation, dot_p - this->length);
}
return {violation, violation <= 0};
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Real Cylinder<DIM>::min_in_direction(const Real * dir, const Real & unit) const {
Real norm = -vector_norm<DIM>(dir);
Real dist = this->in_direction(dir, norm);
return dist/unit;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Real Cylinder<DIM>::max_in_direction(const Real * dir, const Real & unit) const {
Real norm = vector_norm<DIM>(dir);
Real dist = this->in_direction(dir, norm);
return dist/unit;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Real Cylinder<DIM>::size_in_direction(const Real * dir, const Real & unit) const {
Real norm = vector_norm<DIM>(dir);
Real dist = this->in_direction(dir, norm) - this->in_direction(dir, -norm);
return dist/unit;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Cylinder<DIM>::generate_surface_points(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real step_size[DIM],
const Real & repr_const,
const bool periodicity[DIM],
PointContainer<DIM> & points) const
{
Real err2 = 0;
bool periodicity_exists = false;
for (Uint i = 0 ; i < DIM ; ++i) {
err2 += this->axis[i] - static_cast<Real>(periodicity[i]);
periodicity_exists |= periodicity[i];
}
if (periodicity_exists && (err2 > 8*std::numeric_limits<Real>::epsilon())) {
std::stringstream error;
error << "Periodicity (if any) MUST be aligned with the cylinder axis. "
<< "This axis is ";
print_array(this->axis, DIM, "axis", error);
error << std::endl << " and the periodicity is ";
print_array(periodicity, DIM, "periodicity", error);
throw CylinderException(error);
-
}
+
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
inline Real circle_extreme(const Real * centre, const Real * axis_unit,
const Real * dir_norm, const Real & radius) {
+ std::stringstream error;
+ error << "circle_extreme is not implemented for DIM = " << DIM << ".";
+ throw CylinderException(error);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+inline Real circle_extreme<threeD>(const Real * centre, const Real * axis_unit,
+ const Real * dir_norm, const Real & radius) {
+ const Uint DIM = threeD;
Real helper[DIM], r1_unit[DIM], r2_unit[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
helper[(i+1)%DIM] = axis_unit[i];
}
cross_prod(axis_unit, helper, r1_unit);
cross_prod(axis_unit, r1_unit, r2_unit);
- std::cout << "norm of r1 = " << vector_norm(r1_unit) << " (should be 1.)"
- << "norm of r2 = " << vector_norm(r2_unit) << " (should be 1.)"
+ std::cout << "norm of r1 = " << vector_norm<DIM>(r1_unit) << " (should be 1.)"
+ << "norm of r2 = " << vector_norm<DIM>(r2_unit) << " (should be 1.)"
<< std::endl;
- Real cos1 = dot_prod(dir_norm, r1_unit);
- Real cos2 = dot_prod(dir_norm, r2_unit);
+ Real cos1 = dot_prod<DIM>(dir_norm, r1_unit);
+ Real cos2 = dot_prod<DIM>(dir_norm, r2_unit);
Real extreme[DIM];
for (Uint i = 0 ; i < DIM ; ++i) {
extreme[i] = centre[i] + radius * (r1_unit[i] * cos1 + r2_unit[i] * cos2);
}
- return dot_prod(extreme, dir_norm);
+ return dot_prod<DIM>(extreme, dir_norm);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+inline Real circle_extreme<twoD>(const Real * centre, const Real * axis_unit,
+ const Real * dir_norm, const Real & radius) {
+ const Uint DIM = twoD;
+ Real helper[DIM];
+ vector_set<DIM>(helper, centre);
+ vector_incr<DIM>(helper, dir_norm, radius);
+ return dot_prod<DIM>(helper, dir_norm);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Real Cylinder<DIM>::in_direction(const Real * dir, const Real & norm) const {
Real dir_norm[DIM];
vector_set<DIM>(dir_norm, dir);
vector_scale<DIM>(dir_norm, 1./norm);
// check base
- Real dist_base = circle_extreme(this->base, this->axis, dir_norm, this->radius);
+ Real dist_base = circle_extreme<DIM>(this->base.getArray(), this->axis,
+ dir_norm, this->radius);
+ if (DIM == twoD) {
+ return dist_base;
+ }
// check top
Real top[DIM];
- vector_set<DIM>(top, this->base);
+ vector_set<DIM>(top, this->base.getArray());
vector_incr<DIM>(top, this->axis, this->length);
- Real dist_top = circle_extreme(top, this->axis, dir_norm, this->radius);
+ Real dist_top = circle_extreme<DIM>(top, this->axis, dir_norm, this->radius);
if (norm*dist_base > norm*dist_top) {
return dist_base;
}
return dist_top;
}
+
+template class Cylinder<twoD>;
+template class Cylinder<threeD>;
diff --git a/src/geometry/geometry.hh b/src/geometry/geometry.hh
index 336bd02..67027c1 100644
--- a/src/geometry/geometry.hh
+++ b/src/geometry/geometry.hh
@@ -1,181 +1,190 @@
/**
* @file geometry.hh
* @author Till Junge <junge@lsmspc42.epfl.ch>
* @date Wed Apr 4 15:15:10 2012
*
* @brief Geometry motherclass (interface) for CADD mesher
*
* @section LICENSE
*
* <insert lisence here>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __GEOMETRY_HH__
#define __GEOMETRY_HH__
#include "common.hh"
#include "point_container.hh"
#include <sstream>
//#include "cadd_mesh.hh"
template <Uint DIM>
class Mesh;
struct InsideObject {
Real distance;
bool is_inside;
};
/*! \brief Geometry interface for CADD mesher
* You will never use this class.*/
template <Uint DIM>
class Geometry {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/*! Constructor
*\param _name a string identifier for the geometry */
Geometry(std::string _name = ""):name(_name){};
virtual ~Geometry(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// function to print the content of the class
virtual void printself(std::ostream & stream, int indent = 0) const
{
stream << "This is just a interface class and this function should never "
<< "ever be called" <<std::endl;
}
virtual void shift(const PointRef<DIM> & offset) = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
/// returns the string identifier
const std::string & get_name() const
{
return this->name;
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
public:
virtual Geometry * resolveType() const = 0;
/*! test whether a point is inside the domain
* \param point coordinates of the point to test*/
virtual bool is_inside(const Real * point, Real tol) const = 0;
inline bool is_inside(const PointRef<DIM> & point, Real tol) const {
return this->is_inside(&(point.getComponent(0)), tol);}
/*! Returns the minimum value of all constraint violations (negative return
* value means point is inside and return value from the nearest boundary)
* \param point coordintates of the point to test*/
virtual InsideObject from_border(const Real * point, Real tol) const = 0;
inline InsideObject from_border(const PointRef<DIM> & point, Real tol) const {
return this->from_border(&(point.getComponent(0)), tol);}
/*! projects the geometry on a line of direction \a dir and returns the
* minimum of the line segment spanned
* \param dir direction vector (will be scaled to unity by the method
* \param unit unit in which the distance should be expressed, by default unity*/
virtual Real min_in_direction(const Real * dir, const Real & unit = 1) const = 0;
/*! projects the geometry on a line of direction \a dir and returns the
* maximum of the line segment spanned
* \param dir direction vector (will be scaled to unity by the method
* \param unit unit in which the distance should be expressed, by default unity*/
virtual Real max_in_direction(const Real * dir, const Real & unit = 1) const = 0;
/*! projects the geometry on a line of direction \a dir and returns the length
* of the projection
* \param dir direction vector (will be scaled to unity by the method
* \param unit unit in which the distance should be expressed, by default unity*/
virtual Real size_in_direction(const Real * dir, const Real & unit = 1) const = 0;
+ /*! put nodes onto the surface of the geometry
+ * \param auxiliary mesh for refinement interpolation
+ * \param refinement vector of refinements used with the auxiliary mesh
+ * \param step_size array of periodic cell sizes as returned by
+ * lattice::getSpatialConstants()
+ * \param repr_const representative constant of the periodic cell
+ * \param periodicity self-explanatory
+ * \param points container of nodal points for meshing
+ */
virtual void generate_surface_points(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real step_size[DIM],
const Real & repr_const,
const bool periodicity[DIM],
PointContainer<DIM> & points) const = 0;
/*! When the dimensions of the overall geometry are not a power of two times
* the lattice size, the domain can be not filled correctly. This method
* corrects the fucked up side by periodically repeating the ok side.
* \param direction positive values mean that the min side is used to augment
* the max size, negative values the inverse.
* \param dim dimension to fix: 0:x, 1:y, 2:z
*/
virtual void complement_periodically(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
PointContainer<DIM> & points,
int direction,
Uint dim) const = 0;
protected:
std::string name;
};
/// standard output stream operator
template<Uint DIM>
inline std::ostream & operator <<(std::ostream & stream, const Geometry<DIM> & _this)
{
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM, typename subGeom>
Geometry<DIM> * newGeom(const subGeom & geometry) {
return static_cast<Geometry<DIM>*>(new subGeom(geometry));
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
inline bool check_if_ok(const PointRef<DIM> & point,
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real & step_size,
PointContainer<DIM> & points) {
Real ref_sq = square(step_size * auxiliary.interpolate(point,
refinement));
Real step_sq = step_size*step_size;
// if refinement is smaller than 1 => full resolution
if (ref_sq < step_sq) {
return false;
}
// check whether the new point in too close to an existing point
for (Uint point_id = 0; point_id < points.getSize() ; ++point_id) {
if (ref_sq > distance2<DIM>(point, points.getPoint(point_id))) {
return false;
}
}
return true;
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
inline bool add_if_ok(const PointRef<DIM> & point,
const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const Real & step_size,
PointContainer<DIM> & points) {
if (check_if_ok(point, auxiliary, refinement, step_size, points)) {
points.addPoint(point);
return true;
} else {
return false;
}
}
#endif /* __GEOMETRY_HH__ */