mirror of
https://github.com/NGSolve/netgen.git
synced 2024-11-11 16:49:16 +05:00
251 lines
5.2 KiB
C++
251 lines
5.2 KiB
C++
|
#include <mystdlib.h>
|
||
|
|
||
|
#include <myadt.hpp>
|
||
|
#include <gprim.hpp>
|
||
|
|
||
|
namespace netgen
|
||
|
{
|
||
|
|
||
|
ostream & operator<<(ostream & s, const Point2d & p)
|
||
|
{
|
||
|
return s << "(" << p.px << ", " << p.py << ")";
|
||
|
}
|
||
|
|
||
|
ostream & operator<<(ostream & s, const Vec2d & v)
|
||
|
{
|
||
|
return s << "(" << v.vx << ", " << v.vy << ")";
|
||
|
}
|
||
|
|
||
|
#ifdef none
|
||
|
ostream & operator<<(ostream & s, const Line2d & l)
|
||
|
{
|
||
|
return s << l.p1 << "-" << l.p2;
|
||
|
}
|
||
|
|
||
|
ostream & operator<<(ostream & s, const TRIANGLE2D & t)
|
||
|
{
|
||
|
return s << t.p1 << "-" << t.p2 << "-" << t.p3;
|
||
|
}
|
||
|
#endif
|
||
|
|
||
|
|
||
|
double Fastatan2 (double x, double y)
|
||
|
{
|
||
|
if (y > 0)
|
||
|
{
|
||
|
if (x > 0)
|
||
|
return y / (x+y);
|
||
|
else
|
||
|
return 1 - x / (y-x);
|
||
|
}
|
||
|
else if (y < 0)
|
||
|
{
|
||
|
if (x < 0)
|
||
|
return 2 + y / (x+y);
|
||
|
else
|
||
|
return 3 - x / (y-x);
|
||
|
}
|
||
|
else
|
||
|
{
|
||
|
if (x >= 0)
|
||
|
return 0;
|
||
|
else
|
||
|
return 2;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
|
||
|
double Angle (const Vec2d & v)
|
||
|
{
|
||
|
if (v.X() == 0 && v.Y() == 0)
|
||
|
return 0;
|
||
|
|
||
|
double ang = atan2 (v.Y(), v.X());
|
||
|
if (ang < 0) ang+= 2 * M_PI;
|
||
|
return ang;
|
||
|
}
|
||
|
|
||
|
double FastAngle (const Vec2d & v)
|
||
|
{
|
||
|
return Fastatan2 (v.X(), v.Y());
|
||
|
}
|
||
|
|
||
|
double Angle (const Vec2d & v1, const Vec2d & v2)
|
||
|
{
|
||
|
double ang = Angle(v2) - Angle(v1);
|
||
|
if (ang < 0) ang += 2 * M_PI;
|
||
|
return ang;
|
||
|
}
|
||
|
|
||
|
double FastAngle (const Vec2d & v1, const Vec2d & v2)
|
||
|
{
|
||
|
double ang = FastAngle(v2) - FastAngle(v1);
|
||
|
if (ang < 0) ang += 4;
|
||
|
return ang;
|
||
|
}
|
||
|
|
||
|
/*
|
||
|
int CW (const Point2d & p1,const Point2d & p2,const Point2d & p3)
|
||
|
{
|
||
|
return Cross (p2 - p1, p3 - p2) < 0;
|
||
|
}
|
||
|
|
||
|
int CCW (const Point2d & p1,const Point2d & p2,const Point2d & p3)
|
||
|
{
|
||
|
return Cross (p2 - p1, p3 - p2) > 0;
|
||
|
}
|
||
|
*/
|
||
|
|
||
|
double Dist2(const Line2d & g, const Line2d & h )
|
||
|
{
|
||
|
double dd = 0.0, d1,d2,d3,d4;
|
||
|
Point2d cp = CrossPoint(g,h);
|
||
|
|
||
|
if ( Parallel(g,h) || !IsOnLine(g,cp) || !IsOnLine(h,cp) )
|
||
|
{
|
||
|
d1 = Dist2(g.P1(),h.P1());
|
||
|
d2 = Dist2(g.P1(),h.P2());
|
||
|
d3 = Dist2(g.P2(),h.P1());
|
||
|
d4 = Dist2(g.P2(),h.P2());
|
||
|
if (d1<d2) d2 = d1;
|
||
|
if (d3<d4) d4 = d3;
|
||
|
dd = ( d2 < d4 ) ? d2 : d4;
|
||
|
}
|
||
|
return dd;
|
||
|
}
|
||
|
|
||
|
|
||
|
Point2d CrossPoint (const Line2d & l1, const Line2d & l2)
|
||
|
{
|
||
|
double den = Cross (l1.Delta(), l2.Delta());
|
||
|
double num = Cross ( (l2.P1() - l1.P1()), l2.Delta());
|
||
|
|
||
|
if (den == 0)
|
||
|
return l1.P1();
|
||
|
else
|
||
|
return l1.P1() + (num/den) * l1.Delta();
|
||
|
}
|
||
|
|
||
|
|
||
|
int CrossPointBarycentric (const Line2d & l1, const Line2d & l2,
|
||
|
double & lam1, double & lam2)
|
||
|
{
|
||
|
// p = l1.1 + lam1 (l1.2-l1.1) = l2.1 + lam2 (l2.2-l2.1)
|
||
|
double a11 = l1.p2.X() - l1.p1.X();
|
||
|
double a21 = l1.p2.Y() - l1.p1.Y();
|
||
|
double a12 = -(l2.p2.X() - l2.p1.X());
|
||
|
double a22 = -(l2.p2.Y() - l2.p1.Y());
|
||
|
|
||
|
double b1 = l2.p1.X() - l1.p1.X();
|
||
|
double b2 = l2.p1.Y() - l1.p1.Y();
|
||
|
|
||
|
double det = a11*a22 - a12 * a21;
|
||
|
if (det == 0)
|
||
|
return 1;
|
||
|
|
||
|
lam1 = (a22 * b1 - a12 * b2) / det;
|
||
|
lam2 = (a11 * b2 - a21 * b1) / det;
|
||
|
return 0;
|
||
|
}
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
int Parallel (const Line2d & l1, const Line2d & l2, double peps)
|
||
|
{
|
||
|
double p = fabs (Cross (l1.Delta(), l2.Delta()));
|
||
|
// (*mycout) << endl << p << " " << l1.Length() << " " << l2.Length() << endl;
|
||
|
return p <= peps * l1.Length() * l2.Length();
|
||
|
}
|
||
|
|
||
|
int IsOnLine (const Line2d & l, const Point2d & p, double heps)
|
||
|
{
|
||
|
double c1 = (p - l.P1()) * l.Delta();
|
||
|
double c2 = (p - l.P2()) * l.Delta();
|
||
|
double d = fabs (Cross ( (p - l.P1()), l.Delta()));
|
||
|
double len2 = l.Length2();
|
||
|
|
||
|
return c1 >= -heps * len2 && c2 <= heps * len2 && d <= heps * len2;
|
||
|
}
|
||
|
|
||
|
#ifdef none
|
||
|
int IsOnLine (const PLine2d & l, const Point2d & p, double heps)
|
||
|
{
|
||
|
double c1 = (p - l.P1()) * l.Delta();
|
||
|
double c2 = (p - l.P2()) * l.Delta();
|
||
|
double d = fabs (Cross ( (p - l.P1()), l.Delta()));
|
||
|
double len2 = l.Length2();
|
||
|
|
||
|
return c1 >= -heps * len2 && c2 <= heps * len2 && d <= heps * len2;
|
||
|
}
|
||
|
|
||
|
int IsOnLongLine (const Line2d & l, const Point2d & p)
|
||
|
{
|
||
|
double d = fabs (Cross ( (p - l.P1()), l.Delta()));
|
||
|
return d <= EPSGEOM * l.Length();
|
||
|
}
|
||
|
|
||
|
int Hit (const Line2d & l1, const Line2d & l2, double heps)
|
||
|
{
|
||
|
double den = Cross ( (l1.P2() - l1.P1()), (l2.P1() - l2.P2()));
|
||
|
double num1 = Cross ( (l2.P1() - l1.P1()), (l2.P1() - l2.P2()));
|
||
|
double num2 = Cross ( (l1.P2() - l1.P1()), (l2.P1() - l1.P1()));
|
||
|
num1 *= sgn (den);
|
||
|
num2 *= sgn (den);
|
||
|
den = fabs (den);
|
||
|
|
||
|
int ch = (-den * heps <= num1 && num1 <= den * (1 + heps) &&
|
||
|
-den * heps <= num2 && num2 <= den * (1 + heps));
|
||
|
return ch;
|
||
|
}
|
||
|
|
||
|
|
||
|
void Line2d :: GetNormal (Line2d & n) const
|
||
|
{
|
||
|
double ax = P2().X()-P1().X(),
|
||
|
ay = P2().Y()-P1().Y();
|
||
|
Point2d mid(P1().X()+.5*ax, P1().Y()+.5*ay);
|
||
|
|
||
|
n=Line2d(mid,Point2d(mid.X()+ay,mid.Y()-ax)) ;
|
||
|
}
|
||
|
|
||
|
Vec2d Line2d :: NormalDelta () const
|
||
|
{
|
||
|
Line2d tmp;
|
||
|
GetNormal(tmp);
|
||
|
return tmp.Delta();
|
||
|
}
|
||
|
|
||
|
int TRIANGLE2D :: IsOn (const Point2d & p) const
|
||
|
{
|
||
|
return IsOnLine (Line2d (p1, p2), p) ||
|
||
|
IsOnLine (Line2d (p1, p3), p) ||
|
||
|
IsOnLine (Line2d (p2, p3), p);
|
||
|
}
|
||
|
|
||
|
|
||
|
int TRIANGLE2D :: IsIn (const Point2d & p) const
|
||
|
{
|
||
|
return ::CW(p, p1, p2) == ::CW(p, p2, p3) &&
|
||
|
::CW(p, p1, p2) == ::CW(p, p3, p1);
|
||
|
}
|
||
|
|
||
|
|
||
|
|
||
|
int PTRIANGLE2D :: IsOn (const Point2d & p) const
|
||
|
{
|
||
|
return IsOnLine (Line2d (*p1, *p2), p) ||
|
||
|
IsOnLine (Line2d (*p1, *p3), p) ||
|
||
|
IsOnLine (Line2d (*p2, *p3), p);
|
||
|
}
|
||
|
|
||
|
|
||
|
int PTRIANGLE2D :: IsIn (const Point2d & p) const
|
||
|
{
|
||
|
return ::CW(p, *p1, *p2) == ::CW(p, *p2, *p3) &&
|
||
|
::CW(p, *p1, *p2) == ::CW(p, *p3, *p1);
|
||
|
}
|
||
|
|
||
|
#endif
|
||
|
}
|