Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65183684
block.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jun 1, 15:14
Size
12 KB
Mime Type
text/x-c++
Expires
Mon, Jun 3, 15:14 (2 d)
Engine
blob
Format
Raw Data
Handle
18021059
Attached To
rCADDMESH CADD_mesher
block.cc
View Options
#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<twoD>;
template class Block<threeD>;
Event Timeline
Log In to Comment