netgen/libsrc/gprim/geomobjects.hpp

491 lines
10 KiB
C++
Raw Normal View History

2009-01-13 04:40:13 +05:00
#ifndef FILE_OBJECTS
#define FILE_OBJECTS
/* *************************************************************************/
/* File: geomobjects.hpp */
/* Author: Joachim Schoeberl */
/* Date: 20. Jul. 02 */
/* *************************************************************************/
2011-02-28 19:17:25 +05:00
namespace netgen
2009-01-13 04:40:13 +05:00
{
template <int D, typename T = double> class Vec;
template <int D, typename T = double> class Point;
2009-01-13 04:40:13 +05:00
template <int D, typename T>
class Point
2009-01-13 04:40:13 +05:00
{
2011-02-28 19:17:25 +05:00
protected:
T x[D];
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
public:
Point () { ; }
Point (T ax) { for (int i = 0; i < D; i++) x[i] = ax; }
Point (T ax, T ay)
2014-08-20 01:58:36 +06:00
{
// static_assert(D==2, "Point<D> constructor with 2 args called");
x[0] = ax; x[1] = ay;
}
Point (T ax, T ay, T az)
2014-08-20 01:58:36 +06:00
{
// static_assert(D==3, "Point<D> constructor with 3 args called");
x[0] = ax; x[1] = ay; x[2] = az;
}
Point (T ax, T ay, T az, T au)
2011-02-28 19:17:25 +05:00
{ x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;}
2009-01-13 04:40:13 +05:00
2016-10-30 19:01:52 +05:00
template <typename T2>
Point (const Point<D,T2> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2(i); }
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
explicit Point (const Vec<D> & v)
{ for (int i = 0; i < D; i++) x[i] = v(i); }
2009-01-13 04:40:13 +05:00
2016-11-01 15:47:49 +05:00
template <typename T2>
Point & operator= (const Point<D,T2> & p2)
2011-02-28 19:17:25 +05:00
{
2016-11-01 15:47:49 +05:00
for (int i = 0; i < D; i++) x[i] = p2(i);
2011-02-28 19:17:25 +05:00
return *this;
}
2009-01-13 04:40:13 +05:00
Point & operator= (T val)
2011-02-28 19:17:25 +05:00
{
for (int i = 0; i < D; i++) x[i] = val;
return *this;
}
2009-01-13 04:40:13 +05:00
T & operator() (int i) { return x[i]; }
const T & operator() (int i) const { return x[i]; }
2009-01-13 04:40:13 +05:00
T& operator[] (int i) { return x[i]; }
const T& operator[] (int i) const { return x[i]; }
operator const T* () const { return x; }
2018-11-29 22:35:30 +05:00
void DoArchive(Archive& archive)
{
for(int i=0; i<D; i++)
archive & x[i];
}
2011-02-28 19:17:25 +05:00
};
2009-01-13 04:40:13 +05:00
template <int D, typename T>
class Vec
2009-01-13 04:40:13 +05:00
{
2011-02-28 19:17:25 +05:00
protected:
T x[D];
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
public:
Vec () { ; } // for (int i = 0; i < D; i++) x[i] = 0; }
Vec (T ax) { for (int i = 0; i < D; i++) x[i] = ax; }
Vec (T ax, T ay)
2014-08-20 01:58:36 +06:00
{
// static_assert(D==2, "Vec<D> constructor with 2 args called");
x[0] = ax; x[1] = ay;
}
Vec (T ax, T ay, T az)
2014-08-20 01:58:36 +06:00
{
// static_assert(D==3, "Vec<D> constructor with 3 args called");
x[0] = ax; x[1] = ay; x[2] = az;
}
Vec (T ax, T ay, T az, T au)
2011-02-28 19:17:25 +05:00
{ x[0] = ax; x[1] = ay; x[2] = az; x[3] = au; }
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
Vec (const Vec<D> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
2009-01-13 04:40:13 +05:00
2016-10-30 19:01:52 +05:00
explicit Vec (const Point<D,T> & p)
2011-02-28 19:17:25 +05:00
{ for (int i = 0; i < D; i++) x[i] = p(i); }
2009-01-13 04:40:13 +05:00
explicit Vec(const Point<D,T>& p1, const Point<D,T>& p2)
{ for(int i=0; i<D; i++) x[i] = p2(i)-p1(i); }
2009-01-13 04:40:13 +05:00
2016-11-01 15:47:49 +05:00
template <typename T2>
Vec & operator= (const Vec<D,T2> & p2)
2011-02-28 19:17:25 +05:00
{
2016-11-01 15:47:49 +05:00
for (int i = 0; i < D; i++) x[i] = p2(i);
2011-02-28 19:17:25 +05:00
return *this;
}
2009-01-13 04:40:13 +05:00
Vec & operator= (T s)
2011-02-28 19:17:25 +05:00
{
for (int i = 0; i < D; i++) x[i] = s;
return *this;
}
2009-01-13 04:40:13 +05:00
bool operator== (const Vec<D,T> &a) const
{
bool res = true;
for (auto i : Range(D))
res &= (x[i]==a.x[i]);
return res;
}
T & operator() (int i) { return x[i]; }
const T & operator() (int i) const { return x[i]; }
2009-01-13 04:40:13 +05:00
T& operator[] (int i) { return x[i]; }
const T& operator[] (int i) const { return x[i]; }
operator const T* () const { return x; }
2009-01-13 04:40:13 +05:00
2018-11-29 22:35:30 +05:00
void DoArchive(Archive& archive)
{
for(int i=0; i<D; i++)
archive & x[i];
}
T Length () const
2011-02-28 19:17:25 +05:00
{
T l = 0;
2011-02-28 19:17:25 +05:00
for (int i = 0; i < D; i++)
l += x[i] * x[i];
return sqrt (l);
}
2009-01-13 04:40:13 +05:00
T Length2 () const
2011-02-28 19:17:25 +05:00
{
T l = 0;
2011-02-28 19:17:25 +05:00
for (int i = 0; i < D; i++)
l += x[i] * x[i];
return l;
}
2016-10-30 19:01:52 +05:00
Vec & Normalize ()
2011-02-28 19:17:25 +05:00
{
T l = Length();
2016-10-30 19:01:52 +05:00
// if (l != 0)
for (int i = 0; i < D; i++)
x[i] /= (l+1e-40);
2011-02-28 19:17:25 +05:00
return *this;
}
Vec<D> GetNormal () const;
};
2009-01-13 04:40:13 +05:00
2019-09-29 19:22:00 +05:00
template<int D>
inline Vec<D> operator-(const Point<D>& p1, const Point<D>& p2)
{
Vec<D> result;
for(auto i : Range(D))
result[i] = p1[i] - p2[i];
return result;
}
inline double Cross2(const Vec<2>& v1, const Vec<2>& v2)
{
return v1[0] * v2[1] - v1[1] * v2[0];
}
2009-01-13 04:40:13 +05:00
2019-09-29 19:22:00 +05:00
// are points clockwise?
inline bool CW(const Point<2>& p1, const Point<2>& p2,
const Point<2>& p3)
{
return Cross2(p2-p1, p3-p2) < 0;
}
2009-01-13 04:40:13 +05:00
2019-09-29 19:22:00 +05:00
// are points counterclockwise?
inline bool CCW(const Point<2>& p1, const Point<2>& p2,
const Point<2>& p3)
{
return Cross2(p2-p1, p3-p2) > 0;
}
// are strictly points counterclockwise?
inline bool CCW(const Point<2>& p1, const Point<2>& p2,
const Point<2>& p3, double eps)
{
auto v1 = p2-p1;
auto v2 = p3-p2;
return Cross2(v1, v2) > eps*eps*max2(v1.Length2(),
v2.Length2());
}
2009-01-13 04:40:13 +05:00
template <int H, int W=H, typename T = double>
class Mat
2009-01-13 04:40:13 +05:00
{
2011-02-28 19:17:25 +05:00
protected:
T x[H*W];
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
public:
Mat () { ; }
Mat (const Mat & b)
{ for (int i = 0; i < H*W; i++) x[i] = b.x[i]; }
Mat & operator= (T s)
2011-02-28 19:17:25 +05:00
{
for (int i = 0; i < H*W; i++) x[i] = s;
return *this;
}
Mat & operator= (const Mat & b)
{
for (int i = 0; i < H*W; i++) x[i] = b.x[i];
return *this;
}
T & operator() (int i, int j) { return x[i*W+j]; }
const T & operator() (int i, int j) const { return x[i*W+j]; }
T & operator() (int i) { return x[i]; }
const T & operator() (int i) const { return x[i]; }
2011-02-28 19:17:25 +05:00
Vec<H,T> Col (int i) const
2011-02-28 19:17:25 +05:00
{
Vec<H,T> hv;
2011-02-28 19:17:25 +05:00
for (int j = 0; j < H; j++)
hv(j) = x[j*W+i];
return hv;
}
Vec<W,T> Row (int i) const
2011-02-28 19:17:25 +05:00
{
Vec<W,T> hv;
2011-02-28 19:17:25 +05:00
for (int j = 0; j < W; j++)
hv(j) = x[i*W+j];
return hv;
}
void Solve (const Vec<H,T> & rhs, Vec<W,T> & sol) const
2011-02-28 19:17:25 +05:00
{
Mat<W,H,T> inv;
2011-02-28 19:17:25 +05:00
CalcInverse (*this, inv);
sol = inv * rhs;
}
};
2009-01-13 04:40:13 +05:00
2010-09-23 20:07:12 +06:00
2011-02-28 19:17:25 +05:00
template <int D>
class Box
2010-09-23 20:07:12 +06:00
{
2011-02-28 19:17:25 +05:00
protected:
Point<D> pmin, pmax;
public:
Box () { ; }
2010-09-23 20:07:12 +06:00
2011-02-28 19:17:25 +05:00
Box ( const Point<D> & p1)
{
for (int i = 0; i < D; i++)
pmin(i) = pmax(i) = p1(i);
}
2010-09-23 20:07:12 +06:00
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
Box ( const Point<D> & p1, const Point<D> & p2)
{
for (int i = 0; i < D; i++)
{
pmin(i) = min2(p1(i), p2(i));
pmax(i) = max2(p1(i), p2(i));
}
}
2018-07-29 15:03:56 +05:00
Box (const Point<D> & p1, const Point<D> & p2, const Point<D> & p3)
: Box(p1,p2)
{
Add (p3);
}
2011-02-28 19:17:25 +05:00
enum EB_TYPE { EMPTY_BOX = 1 };
Box ( EB_TYPE et )
{
for (int i = 0; i < D; i++)
{
pmin(i) = 1e99;
pmax(i) = -1e99;
}
// pmin = Point<D> (1e99, 1e99, 1e99);
// pmax = Point<D> (-1e99, -1e99, -1e99);
2011-02-28 19:17:25 +05:00
}
const Point<D> & PMin () const { return pmin; }
const Point<D> & PMax () const { return pmax; }
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
void Set (const Point<D> & p)
{ pmin = pmax = p; }
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
void Add (const Point<D> & p)
{
for (int i = 0; i < D; i++)
{
if (p(i) < pmin(i)) pmin(i) = p(i);
2019-02-02 20:24:07 +05:00
/* else */ if (p(i) > pmax(i)) pmax(i) = p(i);
// optimization invalid for empty-box !
2011-02-28 19:17:25 +05:00
}
}
2013-04-03 02:31:02 +06:00
template <typename T1, typename T2>
void Set (const IndirectArray<T1, T2> & points)
{
2019-08-09 21:09:59 +05:00
// Set (points[points.Begin()]);
Set (points[*points.Range().begin()]);
// for (int i = points.Begin()+1; i < points.End(); i++)
for (int i : points.Range().Modify(1,0))
2013-04-03 02:31:02 +06:00
Add (points[i]);
}
template <typename T1, typename T2>
void Add (const IndirectArray<T1, T2> & points)
{
// for (int i = points.Begin(); i < points.End(); i++)
for (int i : points.Range())
2013-04-03 02:31:02 +06:00
Add (points[i]);
}
2011-02-28 19:17:25 +05:00
Point<D> Center () const
{
Point<D> c;
for (int i = 0; i < D; i++)
c(i) = 0.5 * (pmin(i)+pmax(i));
return c;
}
double Diam () const { return Abs (pmax-pmin); }
Point<D> GetPointNr (int nr) const
{
Point<D> p;
for (int i = 0; i < D; i++)
{
p(i) = (nr & 1) ? pmax(i) : pmin(i);
nr >>= 1;
}
return p;
}
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
bool Intersect (const Box<D> & box2) const
{
for (int i = 0; i < D; i++)
if (pmin(i) > box2.pmax(i) ||
pmax(i) < box2.pmin(i)) return 0;
return 1;
}
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
bool IsIn (const Point<D> & p) const
{
for (int i = 0; i < D; i++)
if (p(i) < pmin(i) || p(i) > pmax(i)) return 0;
return 1;
}
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
void Increase (double dist)
{
for (int i = 0; i < D; i++)
{
pmin(i) -= dist;
pmax(i) += dist;
}
}
2018-11-29 22:35:30 +05:00
void DoArchive(Archive& archive)
{ archive & pmin & pmax; }
2009-01-13 04:40:13 +05:00
};
2011-02-28 19:17:25 +05:00
template <int D>
class BoxSphere : public Box<D>
2009-01-13 04:40:13 +05:00
{
2011-02-28 19:17:25 +05:00
protected:
///
Point<D> c;
///
double diam;
///
double inner;
public:
///
BoxSphere () { };
///
BoxSphere (const Box<D> & box)
: Box<D> (box)
{
CalcDiamCenter();
};
///
BoxSphere ( Point<D> apmin, Point<D> apmax )
: Box<D> (apmin, apmax)
{
CalcDiamCenter();
}
///
const Point<D> & Center () const { return c; }
///
double Diam () const { return diam; }
///
double Inner () const { return inner; }
///
void GetSubBox (int nr, BoxSphere & sbox) const
{
for (int i = 0; i < D; i++)
{
if (nr & 1)
{
sbox.pmin(i) = c(i);
sbox.pmax(i) = this->pmax(i);
}
else
{
sbox.pmin(i) = this->pmin(i);
sbox.pmax(i) = c(i);
}
sbox.c(i) = 0.5 * (sbox.pmin(i) + sbox.pmax(i));
nr >>= 1;
}
sbox.diam = 0.5 * diam;
sbox.inner = 0.5 * inner;
}
///
void CalcDiamCenter ()
{
c = Box<D>::Center ();
diam = Dist (this->pmin, this->pmax);
inner = this->pmax(0) - this->pmin(0);
for (int i = 1; i < D; i++)
if (this->pmax(i) - this->pmin(i) < inner)
inner = this->pmax(i) - this->pmin(i);
}
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
};
2009-01-13 04:40:13 +05:00
2018-04-27 11:36:06 +05:00
#ifdef PARALLEL
template <>
inline MPI_Datatype MyGetMPIType<Vec<3, double> > ()
{
static MPI_Datatype MPI_T = 0;
if (!MPI_T)
{
MPI_Type_contiguous ( 3, MPI_DOUBLE, &MPI_T);
MPI_Type_commit ( &MPI_T );
}
return MPI_T;
};
#endif
2009-01-13 04:40:13 +05:00
2011-02-28 19:17:25 +05:00
}
2009-01-13 04:40:13 +05:00
#endif