mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-29 23:30:33 +05:00
429 lines
7.0 KiB
C++
429 lines
7.0 KiB
C++
#ifndef FILE_GEOMOPS
|
|
#define FILE_GEOMOPS
|
|
|
|
/* *************************************************************************/
|
|
/* File: geomops.hpp */
|
|
/* Author: Joachim Schoeberl */
|
|
/* Date: 20. Jul. 02 */
|
|
/* *************************************************************************/
|
|
|
|
|
|
/*
|
|
|
|
Point - Vector operations
|
|
|
|
*/
|
|
|
|
|
|
|
|
|
|
template <class TA, class TB>
|
|
class SumExpr : public VecExpr<SumExpr<TA, TB> >
|
|
{
|
|
const TA a;
|
|
const TB b;
|
|
public:
|
|
SumExpr (const TA aa, const TB ab) : a(aa), b(ab) { ; }
|
|
double operator() (int i) const { return a(i) + b(i); }
|
|
};
|
|
|
|
template <typename TA, typename TB>
|
|
inline SumExpr<TA,TB>
|
|
operator+ (const VecExpr<TA> & a, const VecExpr<TB> & b)
|
|
{
|
|
return SumExpr<TA,TB> (static_cast <const TA&> (a), static_cast <const TB&> (b));
|
|
}
|
|
|
|
/*
|
|
template <int D1, int D2>
|
|
inline SumExpr<const Vec<D1>&, const Vec<D2>&>
|
|
operator+ (const Vec<D1> & a, const Vec<D2> & b)
|
|
{
|
|
return SumExpr<const Vec<D1>&, const Vec<D2>&> (a, b);
|
|
}
|
|
*/
|
|
|
|
|
|
|
|
|
|
|
|
/*
|
|
template <int D>
|
|
inline Vec<D> operator+ (const Vec<D> & a, const Vec<D> & b)
|
|
{
|
|
Vec<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = a(i) + b(i);
|
|
return res;
|
|
}
|
|
*/
|
|
|
|
template <int D>
|
|
inline Point<D> operator+ (const Point<D> & a, const Vec<D> & b)
|
|
{
|
|
Point<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = a(i) + b(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> operator- (const Point<D> & a, const Point<D> & b)
|
|
{
|
|
Vec<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = a(i) - b(i);
|
|
return res;
|
|
}
|
|
|
|
template <int D>
|
|
inline Point<D> operator- (const Point<D> & a, const Vec<D> & b)
|
|
{
|
|
Point<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = a(i) - b(i);
|
|
return res;
|
|
}
|
|
|
|
template <int D>
|
|
inline Vec<D> operator- (const Vec<D> & a, const Vec<D> & b)
|
|
{
|
|
Vec<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = a(i) - b(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> operator* (double s, const Vec<D> & b)
|
|
{
|
|
Vec<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = s * b(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
template <int D>
|
|
inline double operator* (const Vec<D> & a, const Vec<D> & b)
|
|
{
|
|
double sum = 0;
|
|
for (int i = 0; i < D; i++)
|
|
sum += a(i) * b(i);
|
|
return sum;
|
|
}
|
|
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> operator- (const Vec<D> & b)
|
|
{
|
|
Vec<D> res;
|
|
for (int i = 0; i < D; i++)
|
|
res(i) = -b(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
template <int D>
|
|
inline Point<D> & operator+= (Point<D> & a, const Vec<D> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) += b(i);
|
|
return a;
|
|
}
|
|
|
|
|
|
template <int D, typename T>
|
|
inline Point<D> & operator+= (Point<D> & a, const VecExpr<T> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) += b(i);
|
|
return a;
|
|
}
|
|
|
|
template <int D>
|
|
inline Vec<D> & operator+= (Vec<D> & a, const Vec<D> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) += b(i);
|
|
return a;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
template <int D>
|
|
inline Point<D> & operator-= (Point<D> & a, const Vec<D> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) -= b(i);
|
|
return a;
|
|
}
|
|
|
|
template <int D, typename T>
|
|
inline Point<D> & operator-= (Point<D> & a, const VecExpr<T> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) -= b(i);
|
|
return a;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> & operator-= (Vec<D> & a, const Vec<D> & b)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) -= b(i);
|
|
return a;
|
|
}
|
|
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> & operator*= (Vec<D> & a, double s)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) *= s;
|
|
return a;
|
|
}
|
|
|
|
|
|
template <int D>
|
|
inline Vec<D> & operator/= (Vec<D> & a, double s)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
a(i) /= s;
|
|
return a;
|
|
}
|
|
|
|
|
|
|
|
|
|
// Matrix - Vector operations
|
|
|
|
/*
|
|
template <int H, int W>
|
|
inline Vec<H> operator* (const Mat<H,W> & m, const Vec<W> & v)
|
|
{
|
|
Vec<H> res;
|
|
for (int i = 0; i < H; i++)
|
|
{
|
|
res(i) = 0;
|
|
for (int j = 0; j < W; j++)
|
|
res(i) += m(i,j) * v(j);
|
|
}
|
|
return res;
|
|
}
|
|
*/
|
|
|
|
// thanks to VC60 partial template specialization features !!!
|
|
|
|
inline Vec<2> operator* (const Mat<2,2> & m, const Vec<2> & v)
|
|
{
|
|
Vec<2> res;
|
|
for (int i = 0; i < 2; i++)
|
|
{
|
|
res(i) = 0;
|
|
for (int j = 0; j < 2; j++)
|
|
res(i) += m(i,j) * v(j);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
inline Vec<2> operator* (const Mat<2,3> & m, const Vec<3> & v)
|
|
{
|
|
Vec<2> res;
|
|
for (int i = 0; i < 2; i++)
|
|
{
|
|
res(i) = 0;
|
|
for (int j = 0; j < 3; j++)
|
|
res(i) += m(i,j) * v(j);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
inline Vec<3> operator* (const Mat<3,2> & m, const Vec<2> & v)
|
|
{
|
|
Vec<3> res;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
res(i) = 0;
|
|
for (int j = 0; j < 2; j++)
|
|
res(i) += m(i,j) * v(j);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
inline Vec<3> operator* (const Mat<3,3> & m, const Vec<3> & v)
|
|
{
|
|
Vec<3> res;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
res(i) = 0;
|
|
for (int j = 0; j < 3; j++)
|
|
res(i) += m(i,j) * v(j);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*
|
|
template <int H1, int W1, int H2, int W2>
|
|
inline Mat<H1,W2> operator* (const Mat<H1,W1> & a, const Mat<H2,W2> & b)
|
|
{
|
|
Mat<H1,W2> m;
|
|
for (int i = 0; i < H1; i++)
|
|
for (int j = 0; j < W2; j++)
|
|
{
|
|
double sum = 0;
|
|
for (int k = 0; k < W1; k++)
|
|
sum += a(i,k) * b(k, j);
|
|
m(i,j) = sum;
|
|
}
|
|
return m;
|
|
}
|
|
*/
|
|
|
|
inline Mat<2,2> operator* (const Mat<2,2> & a, const Mat<2,2> & b)
|
|
{
|
|
Mat<2,2> m;
|
|
for (int i = 0; i < 2; i++)
|
|
for (int j = 0; j < 2; j++)
|
|
{
|
|
double sum = 0;
|
|
for (int k = 0; k < 2; k++)
|
|
sum += a(i,k) * b(k, j);
|
|
m(i,j) = sum;
|
|
}
|
|
return m;
|
|
}
|
|
|
|
inline Mat<2,2> operator* (const Mat<2,3> & a, const Mat<3,2> & b)
|
|
{
|
|
Mat<2,2> m;
|
|
for (int i = 0; i < 2; i++)
|
|
for (int j = 0; j < 2; j++)
|
|
{
|
|
double sum = 0;
|
|
for (int k = 0; k < 3; k++)
|
|
sum += a(i,k) * b(k, j);
|
|
m(i,j) = sum;
|
|
}
|
|
return m;
|
|
}
|
|
|
|
|
|
inline Mat<3,2> operator* (const Mat<3,2> & a, const Mat<2,2> & b)
|
|
{
|
|
Mat<3,2> m;
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 2; j++)
|
|
{
|
|
double sum = 0;
|
|
for (int k = 0; k < 2; k++)
|
|
sum += a(i,k) * b(k, j);
|
|
m(i,j) = sum;
|
|
}
|
|
return m;
|
|
}
|
|
|
|
inline Mat<3,3> operator* (const Mat<3,3> & a, const Mat<3,3> & b)
|
|
{
|
|
Mat<3,3> m;
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
double sum = 0;
|
|
for (int k = 0; k < 3; k++)
|
|
sum += a(i,k) * b(k, j);
|
|
m(i,j) = sum;
|
|
}
|
|
return m;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
template <int H, int W>
|
|
inline Mat<W,H> Trans (const Mat<H,W> & m)
|
|
{
|
|
Mat<W,H> res;
|
|
for (int i = 0; i < H; i++)
|
|
for (int j = 0; j < W; j++)
|
|
res(j,i) = m(i,j);
|
|
return res;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
template <int D>
|
|
inline ostream & operator<< (ostream & ost, const Vec<D> & a)
|
|
{
|
|
ost << "(";
|
|
for (int i = 0; i < D-1; i++)
|
|
ost << a(i) << ", ";
|
|
ost << a(D-1) << ")";
|
|
return ost;
|
|
}
|
|
|
|
template <int D>
|
|
inline ostream & operator<< (ostream & ost, const Point<D> & a)
|
|
{
|
|
ost << "(";
|
|
for (int i = 0; i < D-1; i++)
|
|
ost << a(i) << ", ";
|
|
ost << a(D-1) << ")";
|
|
return ost;
|
|
}
|
|
|
|
template <int D>
|
|
inline ostream & operator<< (ostream & ost, const Box<D> & b)
|
|
{
|
|
ost << b.PMin() << " - " << b.PMax();
|
|
return ost;
|
|
}
|
|
|
|
template <int H, int W>
|
|
inline ostream & operator<< (ostream & ost, const Mat<H,W> & m)
|
|
{
|
|
ost << "(";
|
|
for (int i = 0; i < H; i++)
|
|
{
|
|
for (int j = 0; j < W; j++)
|
|
ost << m(i,j) << " ";
|
|
ost << endl;
|
|
}
|
|
return ost;
|
|
}
|
|
|
|
|
|
|
|
|
|
#endif
|