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

# ifndef _RHEO_BASIC_POINT_H
# define _RHEO_BASIC_POINT_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
/// 
/// =========================================================================
//
//    basic_point<T> : set of coordinates of type T
//
//    author: Pierre.Saramito@imag.fr
//

/*Class:point
NAME:  @code{point} - vertex of a mesh
@clindex point
@clindex geo
DESCRIPTION:
  Defines geometrical vertex as an array of coordinates.
  This array is also used as a vector of the three dimensional
  physical space.
End:
*/

#include "rheolef/compiler.h"

namespace rheolef { 
//<point:
template <class T>
class basic_point {
    public:

// typedefs:

	typedef size_t size_type;
	typedef T      float_type;

// allocators:

	explicit basic_point () { _x[0] = T();  _x[1] = T();  _x[2] = T(); } 
	
	explicit basic_point (
		const T& x0, 
		const T& x1 = 0, 
		const T& x2 = 0)
			{ _x[0] = x0; _x[1] = x1; _x[2] = x2; }
	
	template <class T1>
	basic_point<T>(const basic_point<T1>& p)
		{ _x[0] = p._x[0]; _x[1] = p._x[1]; _x[2] = p._x[2]; }

	template <class T1>
	basic_point<T>& operator = (const basic_point<T1>& p)
		{ _x[0] = p._x[0]; _x[1] = p._x[1]; _x[2] = p._x[2]; return *this; }

// accessors:

	T& operator[](int i_coord)	        { return _x[i_coord%3]; }
	const T&  operator[](int i_coord) const { return _x[i_coord%3]; }
	T& operator()(int i_coord)	        { return _x[i_coord%3]; }
	const T&  operator()(int i_coord) const { return _x[i_coord%3]; }


// inputs/outputs:

	std::istream& get (std::istream& s, int d = 3)
	{
	    switch (d) {
	    case 1 : _x[1] = _x[2] = 0; return s >> _x[0];
	    case 2 : _x[2] = 0; return s >> _x[0] >> _x[1];
	    default: return s >> _x[0] >> _x[1] >> _x[2];
	    }
	}
	// output
	std::ostream& put (std::ostream& s, int d = 3) const;

// ccomparators: lexicographic order

	template<int d>
	friend bool lexicographically_less (
		const basic_point<T>& a, const basic_point<T>& b) {  
	    for (size_type i = 0; i < d; i++) {
	      if (a[i] < b[i]) return true;
	      if (a[i] > b[i]) return false;
	    }
	    return false; // equality
	}
// algebra:

	friend bool operator == (const basic_point<T>& u, const basic_point<T>& v)
		{ return u[0] == v[0] && u[1] == v[1] && u[2] == v[2]; }

	basic_point<T>& operator+= (const basic_point<T>& v)
		{ _x[0] += v[0]; _x[1] += v[1]; _x[2] += v[2]; return *this; }

	basic_point<T>& operator-= (const basic_point<T>& v)
		{ _x[0] -= v[0]; _x[1] -= v[1]; _x[2] -= v[2]; return *this; }

	basic_point<T>& operator*= (const T& a)
		{ _x[0] *= a; _x[1] *= a; _x[2] *= a; return *this; }

	basic_point<T>& operator/= (const T& a)
		{ _x[0] /= a; _x[1] /= a; _x[2] /= a; return *this; }

	friend basic_point<T> operator+ (const basic_point<T>& u, const basic_point<T>& v)
		{ return basic_point<T> (u[0]+v[0], u[1]+v[1], u[2]+v[2]); }

	friend basic_point<T> operator- (const basic_point<T>& u)
		{ return basic_point<T> (-u[0], -u[1], -u[2]); }

	friend basic_point<T> operator- (const basic_point<T>& u, const basic_point<T>& v)
		{ return basic_point<T> (u[0]-v[0], u[1]-v[1], u[2]-v[2]); }

	template <class T1>
	friend basic_point<T> operator* (const T1& a, const basic_point<T>& u)
		{ return basic_point<T> (a*u[0], a*u[1], a*u[2]); }

	friend basic_point<T> operator* (const basic_point<T>& u, T a)
		{ return basic_point<T> (a*u[0], a*u[1], a*u[2]); }

	friend basic_point<T> operator/ (const basic_point<T>& u, const T& a)
		{ return basic_point<T> (u[0]/a, u[1]/a, u[2]/a); }

	friend basic_point<T> operator/ (const basic_point<T>& u, basic_point<T> v)
		{ return basic_point<T> (u[0]/v[0], u[1]/v[1], u[2]/v[2]); }

	friend basic_point<T> vect (const basic_point<T>& v, const basic_point<T>& w)
		{ return basic_point<T> (
			v[1]*w[2]-v[2]*w[1],
			v[2]*w[0]-v[0]*w[2],
			v[0]*w[1]-v[1]*w[0]); }
// metric:

	// TODO: non-constant metric
	friend T dot (const basic_point<T>& u, const basic_point<T>& v)
		{ return u[0]*v[0]+u[1]*v[1]+u[2]*v[2]; }
	
	friend T norm2 (const basic_point<T>& u)
		{ return dot(u,u); }

	friend T norm (const basic_point<T>& u)
		{ return ::sqrt(norm2(u)); }

	friend T dist2 (const basic_point<T>& x,  const basic_point<T>& y)
		{ return norm2(x-y); }
	
	friend T dist (const basic_point<T>& x,  const basic_point<T>& y)
		{ return norm(x-y); }

	friend T dist_infty (const basic_point<T>& x,  const basic_point<T>& y)
		{ return max(_my_abs(x[0]-y[0]), 
			     max(_my_abs(x[1]-y[1]), _my_abs(x[2]-y[2]))); }

// data:
	T _x[3];
// internal:
protected:
  	static T _my_abs(const T& x) { return (x > T(0)) ? x : -x; }
};
template <class T>
T vect2d (const basic_point<T>& v, const basic_point<T>& w);

template <class T>
T mixt (const basic_point<T>& u, const basic_point<T>& v, const basic_point<T>& w);

// robust(exact) floating point predicates: return value as (0, > 0, < 0)
// formally: orient2d(a,b,x) = vect2d(a-x,b-x)
template <class T>
T orient2d(const basic_point<T>& a, const basic_point<T>& b,
   	const basic_point<T>& x = basic_point<T>());

// formally: orient3d(a,b,c,x) = mixt3d(a-x,b-x,c-x)
template <class T>
T orient3d(const basic_point<T>& a, const basic_point<T>& b,
   	const basic_point<T>& c, const basic_point<T>& x = basic_point<T>());
//>point:

// -------------------------------------------------------------------------------------
// inline'd
// -------------------------------------------------------------------------------------
// input/output
template<class T>
inline
std::ostream&
basic_point<T>::put (std::ostream& s, int d) const
{
    switch (d) {
    case 1 : return s << _x[0];
    case 2 : return s << _x[0] << " " << _x[1];
    default: return s << _x[0] << " " << _x[1] << " " << _x[2];
    }
}
template<class T>
inline
std::istream& 
operator >> (std::istream& s, basic_point<T>& p)
{
    return s >> p[0] >> p[1] >> p[2];
}
template<class T>
inline
std::ostream& 
operator << (std::ostream& s, const basic_point<T>& p)
{
    return s << p[0] << " " << p[1] << " " << p[2];
}


#ifdef _RHEOLEF_HAVE_DOUBLEDOUBLE
template<>
inline
std::ostream& 
basic_point<doubledouble>::put (std::ostream& s, int d) const
{
    switch (d) {
    case 1 : return s << double(_x[0]);
    case 2 : return s << double(_x[0]) << " " << double(_x[1]);
    default: return s << double(_x[0]) << " " << double(_x[1]) << " " << double(_x[2]);
    }
}
#endif // _RHEOLEF_HAVE_DOUBLEDOUBLE
// ----------------------------------------------------------
// point extra: inlined
// ----------------------------------------------------------
#define def_point_function2(f,F)	\
template<class T>			\
inline					\
basic_point<T>				\
f (const basic_point<T>& x)		\
{					\
  basic_point<T> y;			\
  for (size_t i = 0; i < 3; i++)	\
    y[i] = F(x[i]);			\
  return y;				\
}

#define def_point_function(f)	def_point_function2(f,f)

def_point_function(sqr)
def_point_function2(sqrt,::sqrt)
def_point_function(log)
def_point_function(log10)
def_point_function(exp)
#undef def_point_function2
#undef def_point_function

template<class T1, class T2>
inline
basic_point<T1>
operator/ (const T2& a, const basic_point<T1>& x)
{
  basic_point<T1> y;
  for (size_t i = 0; i < 3; i++)
    if (x[i] != 0) y[i] = a/x[i];
  return y;
}
// ===============================
// non-robust predicates
// ===============================
template<class T>
T
orient2d (const basic_point<T>& a, const basic_point<T>& b,
	  const basic_point<T>& x)
{
  error_macro ("robust predicate orient2d: not yet implemented");
  T ax0 = a[0] - x[0];
  T bx0 = b[0] - x[0];
  T ax1 = a[1] - x[1];
  T bx1 = b[1] - x[1];
  return ax0*bx1 - ax1*bx0;
}
template <class T>
T
orient3d (const basic_point<T>& a, const basic_point<T>& b,
	  const basic_point<T>& c, const basic_point<T>& x)
{
  error_macro ("robust predicate orient3d: not yet implemented");
  T ax0 = a[0] - x[0];
  T bx0 = b[0] - x[0];
  T cx0 = c[0] - x[0];
  T ax1 = a[1] - x[1];
  T bx1 = b[1] - x[1];
  T cx1 = c[1] - x[1];
  T ax2 = a[2] - x[2];
  T bx2 = b[2] - x[2];
  T cx2 = c[2] - x[2];
  return ax0 * (bx1 * cx2 - bx2 * cx1)
       + bx0 * (cx1 * ax2 - cx2 * ax1)
       + cx0 * (ax1 * bx2 - ax2 * bx1);
}
// ===============================
// robust predicates : when double
// impletation in predicate.cc
// ===============================
template<>
double
orient2d (const basic_point<double>& a, const basic_point<double>& b,
	  const basic_point<double>& c);
template<>
double
orient3d (const basic_point<double>& a, const basic_point<double>& b,
	  const basic_point<double>& c, const basic_point<double>& d);

// vect2d() and mixt() are deduced from:
template <class T>
inline
T
vect2d (const basic_point<T>& v, const basic_point<T>& w)
{
  return orient2d (v, w);
}
template <class T>
inline
T
mixt (const basic_point<T>& u, const basic_point<T>& v, const basic_point<T>& w)
{
  return orient2d (u, v, w);
}
}// namespace rheolef
#endif /* _RHEO_BASIC_POINT_H */