Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/src/geometry/cylinder.cc
58a104d7bd95dev/incorporate-qhull

Last Change
Raw File

cylinder.cc
View Options

#include "cylinder.hh"
#include <exception>
#include <sstream>
#include <limits>
#include <cmath>
#include <cadd_mesh.hh>
/* -------------------------------------------------------------------------- */
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);
}
this->length = vector_norm<DIM>(&axis_[0]);
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->base[dir] = centre[dir];
this->axis[dir] = axis_[dir]/this->length;
}
this->radius = radius_;
}
/* -------------------------------------------------------------------------- */
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 = 1.;//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;
this->length = vector_norm<DIM>(axis_);
for (Uint dir = 0 ; dir < DIM ; ++dir) {
this->axis[dir] = axis_[dir]/this->length;
}
this->radius = radius;
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
Cylinder<DIM>::Cylinder(const Cylinder & other)
:Geometry<DIM>(other.get_name())
{
vector_set<DIM>(this->base.getArray(), other.base.getArray());
vector_set<DIM>(this->axis, other.axis);
this->length = other.length;
this->radius = other.radius;
}
/* -------------------------------------------------------------------------- */
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>(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>(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<DIM>(relpos.getArray(), axis);
keep_max(violation, -dot_p-tol);
keep_max(violation, dot_p+tol - 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 max_dist = this->in_direction(dir, norm);
Real min_dist = -this->in_direction(dir, -norm);
Real dist = max_dist - min_dist;
return dist/unit;
}
/* -------------------------------------------------------------------------- */
template<Uint DIM>
inline void circle_filler(const Real r1_unit[DIM], const Real r2_unit[DIM],
const Real centre[DIM], const Real & radius,
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) {
Real alpha = 0.;
while (alpha < 2*pi) {
Real periphery[DIM] = {0};
vector_incr<DIM>(periphery, centre);
vector_incr<DIM>(periphery, r1_unit, cos(alpha) * radius);
vector_incr<DIM>(periphery, r2_unit, sin(alpha) * radius);
PointRef<DIM> tmp(periphery);
add_if_ok<DIM>(tmp, auxiliary, refinement, repr_const, points);
// a bit of trig
Real delta_alpha = 2 * acos(sqrt(1-square(repr_const)/
square(radius)));
alpha += delta_alpha;
}
}
/* -------------------------------------------------------------------------- */
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
{
bool consider_periodicity = this->periodicity_check(periodicity);
Real r1_unit[DIM] = { 0.};
r1_unit[0] = 1.;
Real r2_unit[DIM] = { 0.};
r2_unit[1] = 1.;
if (DIM == twoD) {
circle_filler<DIM>(r1_unit, r2_unit, this->base.getArray(), this->radius,
auxiliary, refinement, step_size, repr_const, periodicity,
points);
} else if (DIM == threeD) {
// compute a vector of length R in the plane of the base
Real helper[DIM] = {0.};
for (Uint i = 0 ; i < DIM ; ++i) {
helper[(i+1)%DIM] = this->axis[i];
}
cross_prod(this->axis, helper, r1_unit);
cross_prod(this->axis, r1_unit, r2_unit);
Real top[DIM] = {0.}, bottom[DIM] = {0.};
vector_incr<DIM>(top, this->axis, this->length);
vector_incr<DIM>(top, this->base.getArray());
vector_set<DIM>(bottom, this->base.getArray());
//start by filling the top and bottom of the cylinder
for (Real fill_radius = this->radius; fill_radius > 0;
fill_radius -= repr_const) {
circle_filler<DIM>(r1_unit, r2_unit, bottom, fill_radius,
auxiliary, refinement, step_size, repr_const,
periodicity, points);
if (!consider_periodicity) {
circle_filler<DIM>(r1_unit, r2_unit, top, fill_radius, auxiliary,
refinement, step_size, repr_const, periodicity,
points);
}
fill_radius -= repr_const;
}
//fill the cylindrical surface
for (Real offset = 0 ; offset < this->length/2. ; offset += repr_const) {
//set the offsets for the filler
vector_incr<DIM>(bottom, this->axis, repr_const);
vector_incr<DIM>(top, this->axis, -repr_const);
// fill the ring
circle_filler<DIM>(r1_unit, r2_unit, bottom, this->radius,
auxiliary, refinement, step_size, repr_const,
periodicity, points);
circle_filler<DIM>(r1_unit, r2_unit, top, this->radius,
auxiliary, refinement, step_size, repr_const,
periodicity, points);
}
} else {
std::stringstream error;
error << "Not defined for DIM = " << DIM << ".";
throw CylinderException(error);
}
points.cleanDuplicates(1);
}
/* -------------------------------------------------------------------------- */
template <Uint DIM>
void Cylinder<DIM>::complement_periodically(const Mesh<DIM> & auxiliary,
const std::vector<Real> & refinement,
const PointContainer<DIM> & points,
PointContainer<DIM> & complement,
int direction,
Uint dim) const {
bool periodicity[DIM] = {0};
periodicity[dim] = true;
this->periodicity_check(periodicity);
Real offset[DIM] = {0};
/* positive direction means that bottom points are replicated at the top and
* vice-versa */
offset[dim] = std::copysign(this->length, direction);
Real reference_point[DIM];
vector_set<DIM>(reference_point, this->base.getArray());
if (direction < 0) {
vector_incr<DIM>(reference_point, offset, 1.);
}
vector_incr<DIM>(offset, reference_point);
Real tol = std::max(fabs(reference_point[dim]), 1.);
tol = .1;//2000*tol*std::numeric_limits<Real>::epsilon();
Uint nb_points = points.getSize();
PointRef<DIM> point;
for (Uint i = 0; i < nb_points ; ++i) {
point = points.getConstPoint(i);
if (fabs(point.getComponent(dim) - reference_point[dim]) < tol) {
point[dim] = offset[dim];
add_periodic_complement(auxiliary, refinement, point, points, complement);
}
}
}
/* -------------------------------------------------------------------------- */
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);
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<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<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.getArray());
vector_incr<DIM>(top, this->axis, this->length);
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 <Uint DIM>
bool Cylinder<DIM>::periodicity_check(const bool periodicity[DIM]) 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);
}
return periodicity_exists;
}
template class Cylinder<twoD>;
template class Cylinder<threeD>;
c4science · Help