librheolef-dev 5.93-2 / usr / include / rheolef / hazel.h

#ifndef _RHEO_HAZEL_H
#define _RHEO_HAZEL_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
// 
// Subdivision tree on a mesh with grid stump 
// and octree-like embranching. 
//
// author: jocelyn.etienne@imag.fr
//
// date: 14 Feb 2002
// 
#include "rheolef/meshpoint.h"
#include "rheolef/Vector.h"
namespace rheolef { 

#ifdef _RHEOLEF_HAVE_SLIST
# include <slist>
#elif defined(_RHEOLEF_HAVE_SLIST_H)
# include <slist.h>
#else
# include <list>
# ifndef slist
#  define slist list
# endif
#endif

//! hazel tree node basic container 
/*! Does not need to provide neighbourhood information
 */
//<hazel_stump
class hazel_stump
 {
public:
 
	typedef meshpoint::size_type size_type;

// allocators/deallocators
	
	hazel_stump (int level, int* i_stump, int* i_bough);
	~hazel_stump ();

// modifiers
	
	void
	add_interior_element (size_type element_index);

	void
	add_boundary_element (size_type element_index);

	size_type
	pop_interior_element ();

	size_type
	pop_boundary_element ();

// accessors

	//! Depth in the hazel tree, beginning 0
	int
	level () const;

	//! Indices of stump forebear
	void
	stump_coordinates (int* i_stump) const;

	//! Indices in a virtual grid division by 2^level x 2^level x 2^level
	//! of stump forebear : allows calculation of physical coordinates 
	void
	bough_coordinates (int* i_bough) const;

	size_type
	n_boundary_elements () const;

	size_type
	n_interior_elements () const;

	std::slist<size_type>::const_iterator
	begin_interior () const;

	std::slist<size_type>::const_iterator
	end_interior () const;

	std::slist<size_type>::const_iterator
	begin_boundary () const;

	std::slist<size_type>::const_iterator
	end_boundary () const;

	friend std::ostream&
	operator << (std::ostream& s, const hazel_stump& h);

private:
	int i_stump_[3];
	int i_bough_[3];
	int	   level_;

	size_type n_interior_elements_;
	size_type n_boundary_elements_;
	//! List of elements in the leaf
	/*! List of elements intersecting the slab defined by the
	 *  tree leaf. NULL if not a leaf.
	 */
	std::slist<size_type> interior_element_list;
	std::slist<size_type> boundary_element_list;
 };
//>hazel_stump

inline
int
hazel_stump::level () const
 {
	return level_;
 }

inline
void
hazel_stump::stump_coordinates (int* i_stump) const
 {
	for (int i=0;i<3;i++) i_stump[i]=i_stump_[i];
 }

inline
void
hazel_stump::bough_coordinates (int* i_bough) const
 {
	for (int i=0;i<3;i++) i_bough[i]=i_bough_[i];
 }


inline
hazel_stump::size_type
hazel_stump::n_interior_elements () const
 {
	return n_interior_elements_;
 }


inline
hazel_stump::size_type
hazel_stump::n_boundary_elements () const
 {
	return n_boundary_elements_;
 }


inline
std::slist<hazel_stump::size_type>::const_iterator
hazel_stump::begin_boundary () const
 {
 	return boundary_element_list.begin();
 };

inline
std::slist<hazel_stump::size_type>::const_iterator
hazel_stump::end_boundary () const
 {
 	return boundary_element_list.end();
 };

inline
std::slist<hazel_stump::size_type>::const_iterator
hazel_stump::begin_interior () const
 {
 	return interior_element_list.begin();
 };

inline
std::slist<hazel_stump::size_type>::const_iterator
hazel_stump::end_interior () const
 {
 	return interior_element_list.end();
 };

inline
std::ostream&
operator << (std::ostream& s, const hazel_stump& h) 
 {
 	typedef hazel_stump::size_type size_type;
	std::slist<size_type>::const_iterator e =h.boundary_element_list.begin();
	std::slist<size_type>::const_iterator last =h.boundary_element_list.end();
	while (e!=last)
	 {
		s << (*e) << " ";
		e++;
	 };
	return s << "\n";
 };

class hazel_node;

//! hazel tree node (other than a stump) container
/*! These contain additionnal data on neighbouring nodes
 *  (in a given direction, node B is a neighbour of A iff 
 *  B is the deepest node in the tree which A has one side 
 *  fully in common with) 
 */
//<hazel_bough
class hazel_bough:public hazel_stump
 {
public:
	
// allocators/dellocators
	hazel_bough (int l, int * is, int * ib) : hazel_stump (l,is,ib)
		{};
// modifiers

	//! Roughly set the neighbours to the guess raw_neighbours
	/*! raw_neighbours are supposed exact or may be the father's 
	 *  neighbours in directions marked by squirrel_leap.
	 */
	void
	set_neighbours (hazel_node* raw_neighbours[6], 
					const int squirrel_leap[3]);
	
	//! Refines the guess to exact neighbours
	void
	identify_neighbours ();

// accessors

	hazel_node**
	get_neighbours (); 

private:
	/*! Info for acceleration of neighbour-finding */
	int 	squirrel_leap_[3];

	//! List of neighbours
	/*! List of deepest neighbouring nodes among shallower or as deep
	 *  hazel nodes in every direction.
	 */
	hazel_node*	neighbour_[6];
 };
//>hazel_bough;

//! Node of a hazel tree
class hazel_node : public hazel_bough
 {
public:
	hazel_node(int l, int * is, int * ib) : hazel_bough(l, is, ib)
		{};
private:
	hazel_bough** branch;
 };

//! Subdivision tree on a mesh with grid stump 
//! and octree-like embranching. 
//<hazel
class hazel : private Vector<hazel_node*>
 {
public:

// allocators/deallocators
	
	/*! Creates a hazel with a maximal number of elements per leaf as far 
	 * as max_depth allows (if set). An upper bound for the number of stumps 
	 * can be set, which will be approached as close as possible.
	 *
	 * tolerance<0 for using default value (non-zero)
	 */
	hazel 	();
	
	~hazel	();

// initialization

	void
	grow
	 (
		std::vector<geo_element>::const_iterator begin_element,
		std::vector<geo_element>::const_iterator end_element,
		const std::vector<point>::const_iterator begin_node,
		const std::vector<point>::const_iterator end_node,
		std::vector<geo_element>::const_iterator begin_boundary,
		std::vector<geo_element>::const_iterator end_boundary,
		const int mesh_size,
		const int space_dimension,
		const point& x_min, const point& x_max,
		Float tolerance,
		const size_type elements_per_leaf, 
		const size_type max_depth=std::numeric_limits<size_type>::max()
	 );

// definitions

	typedef enum
	 {
	 	left 	= 0,
		right 	= 1,
		down	= 2,
		up		= 3,
		back	= 4,
		forth	= 5
	 } 	direction;
	 
// accessors

	int
	avg_list_size ();
	
	/*! True iff x is in omega
	 *  - In: x[omega.dimension] 
	 *  - Out: element containing x, hazel leaf
	 *  	containing element.
	 */
	
	bool
	localize	(const point& x, geo_element::size_type& element) const;

	/*! - In: x[omega.dimension]
	 *  - Out: y, wich is either :
	 *			- x if x is in omega
	 *			- the projection of x on the nearest segment of omega.
	 *		element containing y, hazel leaf containing element.
	 */
	
	void
	localize_nearest	(const point& x, point& y, geo_element::size_type& element) const;

	/*!
	 *	In: x0 point of omega, v vector.
	 *
	 *	Out:
	 *		- x = x0 + t v with t<=1 such that ]x0 x[ \in omega 
	 *		with  |x0 x| to the maximum.
	 *		- element contains x.
	 */
	bool
	trace (const point& x0, const point& v,
			point& x, Float& t, size_type& element) const;

private:
	// methods :

	bool
	localize_in_leaf (const point& x, hazel_node& leaf, 
		geo_element::size_type& element) const;

	bool in_bbox_element(const geo_element&, const point&) const;
	bool try_hard_in_leaf(const point&, hazel_node&, geo_element::size_type&) const;
	void list_leaf (std::ofstream& os, hazel_node& leaf) const;

	bool
	localize_nearest_in_leaf (const point& x, point& y, const hazel_node& leaf) const;

	bool
	in_reference_element (const geo_element::enum_type t, const point& x) const;

	bool
	in_interior_element (const geo_element& e, const point& x, bool& not_too_far) const;

	bool
	in_interior_element (const geo_element& e, const point& x) const;

	/*! When true, segment [x0 x0+v] hits the element e of dimension d-1
	 * at x = x0 + t v. (out: x and t) 
	 * inward set iff vector v hits from the outside.
	 */
	
	bool
	hit_boundary_element (const geo_element& e,
						  const point& x0, const point& v, 
						  point& x, Float& ti, bool& inward) const;
	
	/*! When true,  segment [x0x1] hits the box xmin,xmax and gives the
	 * parameter t for the point of [x0x1] where the segment enters the
	 * box. "side" and "dir" give the normal of the face of the box 
	 * the segments enters thro. (side=1 and dir=0 mean n = -e1, i.e. 
	 * y-backward). 
	 *
	 * If point x1 belongs to box, then t=1, side=dir=-1.
	 */
	bool
	hit_boundingbox (const point& x0, const point& v, 
					const point& xmin, const point& xmax, 
					Float& t, int& side, int& dir) const;
	
	bool
	in_boundingbox (const point& x, const point& xmin, const point& xmax) const;

	//! Returns front bottom left point coordinates of a node
	point
	box_min	(hazel_node* leaf) const;
	//! Returns back top right point coordinates of a node
	point
	box_max	(hazel_node* leaf) const;

	//! Accesses to the stump index in Vector coming next in dir direction 
	size_type 
	get_neighbour	(const size_type stump, const direction dir) const;

	//! Returns the (linear) index number for given integer coordinates in
	//! 3d grid.
	size_type 
	index	(int* pos) const;
	size_type 
	index	(int posx, int posy, int posz) const;

	int dimension;
	point xmin, xmax;
	std::vector<geo_element>::const_iterator begin_element_;
	std::vector<geo_element>::const_iterator end_element_;
	std::vector<point>::const_iterator begin_node_;
	std::vector<geo_element>::const_iterator begin_boundary_;
	std::vector<geo_element>::const_iterator end_boundary_;

	int n_stump[3];
	point h_stump;
	//! Tolerance for inclusion tests
	Float	tol;

	bool grown;

 };
//>hazel

inline
meshpoint::size_type
hazel::index	(int posx, int posy, int posz) const
 {
	return posx +n_stump[0] *(posy +n_stump[2]*posz);
 };

inline
meshpoint::size_type
hazel::index	(int* pos) const
 {
	return pos[0] +n_stump[0] *(pos[1] +n_stump[2]*pos[2]);
 };

inline
point
hazel::box_min	(hazel_node* leaf) const
 {
	point x;
	int i_stump[3];
	int i_bough[3];
	leaf->stump_coordinates(i_stump);
	leaf->bough_coordinates(i_bough);

	for (int i=0; i<3; i++) x[i] =xmin[i] +i_stump[i]*h_stump[i]
										  +i_bough[i]*h_stump[i]/(1<<leaf->level());
	return x;
 }

inline
point
hazel::box_max	(hazel_node* leaf) const
 {
 	return box_min(leaf)+h_stump/(1<<leaf->level());
 }
inline
bool
hazel::in_boundingbox (const point& x, const point& xmin, const point& xmax) const
{
 	for (int i=0; i<dimension; i++)
	 {
	 	if (x[i]<xmin[i]) return false;
		if (x[i]>xmax[i]) return false;
	 }
	return true;
}
}// namespace rheolef
#endif // _RHEO_HAZEL_H