ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-14 10:08:32 +05:00
Code Issues Packages Projects Releases Wiki Activity

header files

Browse Source
This commit is contained in:
Joachim Schoeberl 2011-02-28 14:17:25 +00:00
parent 9a043aae26
commit 807d091d9e
19 changed files with 2049 additions and 2046 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
libsrc/geom2d/Makefile.am
View File

@ -7,7 +7,6 @@ lib_LTLIBRARIES = libgeom2d.la libgeom2dvis.la
libgeom2d_la_SOURCES = genmesh2d.cpp geom2dmesh.cpp geometry2d.cpp libgeom2d_la_SOURCES = genmesh2d.cpp geom2dmesh.cpp geometry2d.cpp
libgeom2d_la_LIBADD = $(top_builddir)/libsrc/meshing/libmesh.la libgeom2d_la_LIBADD = $(top_builddir)/libsrc/meshing/libmesh.la
# $(top_builddir)/libsrc/gprim/libgprim.la
libgeom2dvis_la_SOURCES = geom2dpkg.cpp vsgeom2d.cpp libgeom2dvis_la_SOURCES = geom2dpkg.cpp vsgeom2d.cpp
libgeom2dvis_la_LIBADD = libgeom2d.la libgeom2dvis_la_LIBADD = libgeom2d.la

6
libsrc/geom2d/genmesh2d.cpp
View File

@ -1,8 +1,6 @@
#include <mystdlib.h>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
namespace netgen namespace netgen
{ {
@ -10,10 +8,6 @@ namespace netgen
void CalcPartition (double l, double h, double h1, double h2, void CalcPartition (double l, double h, double h1, double h2,
double hcurve, double elto0, Array<double> & points); double hcurve, double elto0, Array<double> & points);

2
libsrc/geom2d/geom2dmesh.cpp
View File

@ -1,6 +1,4 @@
#include <mystdlib.h>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
namespace netgen namespace netgen

9
libsrc/geom2d/geom2dpkg.cpp
View File

@ -1,15 +1,8 @@
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <incvis.hpp>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
#include <visual.hpp> #include <visual.hpp>
#include "vsgeom2d.hpp"
#include "vsgeom2d.hpp"
// extern "C" int Ng_CSG_Init (Tcl_Interp * interp); // extern "C" int Ng_CSG_Init (Tcl_Interp * interp);

2
libsrc/geom2d/geometry2d.cpp
View File

@ -4,8 +4,6 @@
*/ */
#include <mystdlib.h>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>

3
libsrc/geom2d/geometry2d.hpp
View File

@ -12,10 +12,9 @@
// #include "../gprim/spline.hpp" // #include "../gprim/spline.hpp"
#include "../gprim/splinegeometry.hpp" // #include "../gprim/splinegeometry.hpp"
#include "geom2dmesh.hpp" #include "geom2dmesh.hpp"
namespace netgen namespace netgen
{ {

5
libsrc/geom2d/vsgeom2d.cpp
View File

@ -1,9 +1,4 @@
#include <mystdlib.h>
#include "incvis.hpp"
#include <myadt.hpp>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
#include <visual.hpp> #include <visual.hpp>

9
libsrc/gprim/adtree.hpp
View File

@ -9,13 +9,15 @@
/* *************************************************************************/ /* *************************************************************************/
namespace netgen
{
/** /**
Alternating Digital Tree Alternating Digital Tree
*/ */
#include "../include/mystdlib.h" // #include "../include/mystdlib.h"
#include "../include/myadt.hpp" // #include "../include/myadt.hpp"
class ADTreeNode class ADTreeNode
{ {
@ -478,4 +480,7 @@ public:
const ADTree6 & Tree() const { return *tree; }; const ADTree6 & Tree() const { return *tree; };
}; };
}
#endif #endif

271
libsrc/gprim/geom2d.hpp
View File

@ -7,74 +7,75 @@
/* Date: 5. Aug. 95 */ /* Date: 5. Aug. 95 */
/* *************************************************************************/ /* *************************************************************************/
namespace netgen
{
/* Geometric Algorithms */
/* Geometric Algorithms */
#define EPSGEOM 1E-5 #define EPSGEOM 1E-5
// extern void MyError (const char * ch); // extern void MyError (const char * ch);
class Point2d; class Point2d;
class Vec2d; class Vec2d;
class LINE2D; class LINE2D;
class Line2d; class Line2d;
class PLine2d; class PLine2d;
class TRIANGLE2D; class TRIANGLE2D;
class PTRIANGLE2D; class PTRIANGLE2D;
inline Vec2d operator- (const Point2d & p1, const Point2d & p2); inline Vec2d operator- (const Point2d & p1, const Point2d & p2);
inline Point2d operator- (const Point2d & p1, const Vec2d & v); inline Point2d operator- (const Point2d & p1, const Vec2d & v);
inline Point2d operator+ (const Point2d & p1, const Vec2d & v); inline Point2d operator+ (const Point2d & p1, const Vec2d & v);
inline Point2d Center (const Point2d & p1, const Point2d & p2); inline Point2d Center (const Point2d & p1, const Point2d & p2);
inline void PpSmV (const Point2d & p1, double s, const Vec2d & v, Point2d & p2); inline void PpSmV (const Point2d & p1, double s, const Vec2d & v, Point2d & p2);
inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v); inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v);
ostream & operator<<(ostream & s, const Point2d & p); ostream & operator<<(ostream & s, const Point2d & p);
inline Vec2d operator- (const Point2d & p1, const Point2d & p2); inline Vec2d operator- (const Point2d & p1, const Point2d & p2);
inline Point2d operator- (const Point2d & p1, const Vec2d & v); inline Point2d operator- (const Point2d & p1, const Vec2d & v);
inline Point2d operator+ (const Point2d & p1, const Vec2d & v); inline Point2d operator+ (const Point2d & p1, const Vec2d & v);
inline Vec2d operator- (const Vec2d & p1, const Vec2d & v); inline Vec2d operator- (const Vec2d & p1, const Vec2d & v);
inline Vec2d operator+ (const Vec2d & p1, const Vec2d & v); inline Vec2d operator+ (const Vec2d & p1, const Vec2d & v);
inline Vec2d operator* (double scal, const Vec2d & v); inline Vec2d operator* (double scal, const Vec2d & v);
double Angle (const Vec2d & v); double Angle (const Vec2d & v);
double FastAngle (const Vec2d & v); double FastAngle (const Vec2d & v);
double Angle (const Vec2d & v1, const Vec2d & v2); double Angle (const Vec2d & v1, const Vec2d & v2);
double FastAngle (const Vec2d & v1, const Vec2d & v2); double FastAngle (const Vec2d & v1, const Vec2d & v2);
ostream & operator<<(ostream & s, const Vec2d & v); ostream & operator<<(ostream & s, const Vec2d & v);
double Dist2(const Line2d & g, const Line2d & h ); // GH double Dist2(const Line2d & g, const Line2d & h ); // GH
int Near (const Point2d & p1, const Point2d & p2, const double eps); int Near (const Point2d & p1, const Point2d & p2, const double eps);
int Parallel (const Line2d & l1, const Line2d & l2, double peps = EPSGEOM); int Parallel (const Line2d & l1, const Line2d & l2, double peps = EPSGEOM);
int IsOnLine (const Line2d & l, const Point2d & p, double heps = EPSGEOM); int IsOnLine (const Line2d & l, const Point2d & p, double heps = EPSGEOM);
int IsOnLongLine (const Line2d & l, const Point2d & p); int IsOnLongLine (const Line2d & l, const Point2d & p);
int Hit (const Line2d & l1, const Line2d & l2, double heps = EPSGEOM); int Hit (const Line2d & l1, const Line2d & l2, double heps = EPSGEOM);
ostream & operator<<(ostream & s, const Line2d & l); ostream & operator<<(ostream & s, const Line2d & l);
Point2d CrossPoint (const PLine2d & l1, const PLine2d & l2); Point2d CrossPoint (const PLine2d & l1, const PLine2d & l2);
Point2d CrossPoint (const Line2d & l1, const Line2d & l2); Point2d CrossPoint (const Line2d & l1, const Line2d & l2);
int Parallel (const PLine2d & l1, const PLine2d & l2, double peps = EPSGEOM); int Parallel (const PLine2d & l1, const PLine2d & l2, double peps = EPSGEOM);
int IsOnLine (const PLine2d & l, const Point2d & p, double heps = EPSGEOM); int IsOnLine (const PLine2d & l, const Point2d & p, double heps = EPSGEOM);
int IsOnLongLine (const PLine2d & l, const Point2d & p); int IsOnLongLine (const PLine2d & l, const Point2d & p);
int Hit (const PLine2d & l1, const Line2d & l2, double heps = EPSGEOM); int Hit (const PLine2d & l1, const Line2d & l2, double heps = EPSGEOM);
ostream & operator<<(ostream & s, const Line2d & l); ostream & operator<<(ostream & s, const Line2d & l);
ostream & operator<<(ostream & s, const TRIANGLE2D & t); ostream & operator<<(ostream & s, const TRIANGLE2D & t);
ostream & operator<<(ostream & s, const PTRIANGLE2D & t); ostream & operator<<(ostream & s, const PTRIANGLE2D & t);
double Dist2 (const Point2d & p1, const Point2d & p2); double Dist2 (const Point2d & p1, const Point2d & p2);
/// ///
class Point2d class Point2d
{ {
/// ///
friend class Vec2d; friend class Vec2d;
protected: protected:
/// ///
double px, py; double px, py;
public: public:
/// ///
Point2d() { /* px = py = 0; */ } Point2d() { /* px = py = 0; */ }
/// ///
@ -192,28 +193,28 @@ public:
/// ///
friend ostream & operator<<(ostream & s, const Point2d & p); friend ostream & operator<<(ostream & s, const Point2d & p);
}; };
inline int Near (const Point2d & p1, const Point2d & p2, inline int Near (const Point2d & p1, const Point2d & p2,
const double eps = 1e-4 ) const double eps = 1e-4 )
{
return Dist2(p1,p2) <= eps*eps;
}
///
class Vec2d
{ {
protected: return Dist2(p1,p2) <= eps*eps;
}
///
class Vec2d
{
protected:
/// ///
double vx, vy; double vx, vy;
public: public:
/// ///
Vec2d() { /* vx = vy = 0; */ } Vec2d() { /* vx = vy = 0; */ }
/// ///
@ -286,7 +287,7 @@ public:
/// ///
friend inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v); friend inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v);
/// Angle in [0,2*PI) /// Angle in [0,2*PI)
/// ///
friend double Angle (const Vec2d & v); friend double Angle (const Vec2d & v);
@ -303,14 +304,14 @@ public:
/// ///
class Line2d class Line2d
{ {
protected: protected:
/// ///
Point2d p1, p2; Point2d p1, p2;
public: public:
/// ///
Line2d() : p1(), p2() { }; Line2d() : p1(), p2() { };
/// ///
@ -375,14 +376,14 @@ public:
#ifdef NONE #ifdef NONE
/// ///
class PLine2d class PLine2d
{ {
protected: protected:
/// ///
Point2d const * p1, *p2; Point2d const * p1, *p2;
public: public:
/// ///
PLine2d() { }; PLine2d() { };
/// ///
@ -440,8 +441,8 @@ public:
/// ///
class ILINE class ILINE
{ {
/// ///
INDEX i[2]; INDEX i[2];
@ -484,14 +485,14 @@ class ILINE
/// ///
class TRIANGLE2D class TRIANGLE2D
{ {
private: private:
/// ///
Point2d p1, p2, p3; Point2d p1, p2, p3;
public: public:
/// ///
TRIANGLE2D() { }; TRIANGLE2D() { };
/// ///
@ -545,14 +546,14 @@ public:
}; };
/// ///
class PTRIANGLE2D class PTRIANGLE2D
{ {
private: private:
/// ///
Point2d const *p1, *p2, *p3; Point2d const *p1, *p2, *p3;
public: public:
/// ///
PTRIANGLE2D() { }; PTRIANGLE2D() { };
/// ///
@ -609,12 +610,12 @@ public:
class Polygon2d class Polygon2d
{ {
protected: protected:
Array<Point2d> points; Array<Point2d> points;
public: public:
Polygon2d (); Polygon2d ();
~Polygon2d (); ~Polygon2d ();
@ -637,32 +638,32 @@ public:
int IsStarPoint (const Point2d & p) const; int IsStarPoint (const Point2d & p) const;
Point2d Center() const; Point2d Center() const;
Point2d EqualAreaPoint () const; Point2d EqualAreaPoint () const;
private: private:
double HArea () const; double HArea () const;
}; };
/** Cheap approximation to atan2. /** Cheap approximation to atan2.
A monotone function of atan2(x,y) is computed. A monotone function of atan2(x,y) is computed.
*/ */
extern double Fastatan2 (double x, double y); extern double Fastatan2 (double x, double y);
inline Vec2d & Vec2d :: operator+= (const Vec2d & v2) inline Vec2d & Vec2d :: operator+= (const Vec2d & v2)
{ {
vx += v2.vx; vx += v2.vx;
vy += v2.vy; vy += v2.vy;
return *this; return *this;
} }
inline Vec2d & Vec2d :: operator-= (const Vec2d & v2) inline Vec2d & Vec2d :: operator-= (const Vec2d & v2)
{ {
vx -= v2.vx; vx -= v2.vx;
vy -= v2.vy; vy -= v2.vy;
return *this; return *this;
} }
inline Vec2d & Vec2d :: operator*= (double s) inline Vec2d & Vec2d :: operator*= (double s)
{ {
vx *= s; vx *= s;
vy *= s; vy *= s;
@ -670,8 +671,8 @@ inline Vec2d & Vec2d :: operator*= (double s)
} }
inline Vec2d & Vec2d :: operator/= (double s) inline Vec2d & Vec2d :: operator/= (double s)
{ {
if (s != 0) if (s != 0)
{ {
vx /= s; vx /= s;
@ -682,53 +683,53 @@ inline Vec2d & Vec2d :: operator/= (double s)
MyError ("Vec2d::operator /=: Division by zero"); MyError ("Vec2d::operator /=: Division by zero");
} }
return *this; return *this;
} }
inline Vec2d operator- (const Point2d & p1, const Point2d & p2) inline Vec2d operator- (const Point2d & p1, const Point2d & p2)
{ {
return Vec2d (p1.X() - p2.X(), p1.Y() - p2.Y()); return Vec2d (p1.X() - p2.X(), p1.Y() - p2.Y());
} }
inline Point2d operator- (const Point2d & p1, const Vec2d & v) inline Point2d operator- (const Point2d & p1, const Vec2d & v)
{ {
return Point2d (p1.X() - v.X(), p1.Y() - v.Y()); return Point2d (p1.X() - v.X(), p1.Y() - v.Y());
} }
inline Point2d operator+ (const Point2d & p1, const Vec2d & v) inline Point2d operator+ (const Point2d & p1, const Vec2d & v)
{ {
return Point2d (p1.X() + v.X(), p1.Y() + v.Y()); return Point2d (p1.X() + v.X(), p1.Y() + v.Y());
} }
inline Point2d Center (const Point2d & p1, const Point2d & p2) inline Point2d Center (const Point2d & p1, const Point2d & p2)
{ {
return Point2d ((p1.X() + p2.X()) / 2, (p1.Y() + p2.Y()) / 2); return Point2d ((p1.X() + p2.X()) / 2, (p1.Y() + p2.Y()) / 2);
} }
inline Vec2d operator- (const Vec2d & v1, const Vec2d & v2) inline Vec2d operator- (const Vec2d & v1, const Vec2d & v2)
{ {
return Vec2d (v1.X() - v2.X(), v1.Y() - v2.Y()); return Vec2d (v1.X() - v2.X(), v1.Y() - v2.Y());
} }
inline Vec2d operator+ (const Vec2d & v1, const Vec2d & v2) inline Vec2d operator+ (const Vec2d & v1, const Vec2d & v2)
{ {
return Vec2d (v1.X() + v2.X(), v1.Y() + v2.Y()); return Vec2d (v1.X() + v2.X(), v1.Y() + v2.Y());
} }
inline Vec2d operator* (double scal, const Vec2d & v) inline Vec2d operator* (double scal, const Vec2d & v)
{ {
return Vec2d (scal * v.X(), scal * v.Y()); return Vec2d (scal * v.X(), scal * v.Y());
} }
inline void PpSmV (const Point2d & p1, double s, inline void PpSmV (const Point2d & p1, double s,
const Vec2d & v, Point2d & p2) const Vec2d & v, Point2d & p2)
{ {
p2.X() = p1.X() + s * v.X(); p2.X() = p1.X() + s * v.X();
@ -736,7 +737,7 @@ inline void PpSmV (const Point2d & p1, double s,
} }
inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v) inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v)
{ {
v.X() = p1.X() - p2.X(); v.X() = p1.X() - p2.X();
v.Y() = p1.Y() - p2.Y(); v.Y() = p1.Y() - p2.Y();
@ -747,19 +748,19 @@ inline void PmP (const Point2d & p1, const Point2d & p2, Vec2d & v)
#ifdef none #ifdef none
inline int TRIANGLE2D :: Regular() const inline int TRIANGLE2D :: Regular() const
{ {
return fabs(Cross ( p2 - p1, p3 - p2)) > EPSGEOM; return fabs(Cross ( p2 - p1, p3 - p2)) > EPSGEOM;
} }
inline int TRIANGLE2D :: CW () const inline int TRIANGLE2D :: CW () const
{ {
return Cross ( p2 - p1, p3 - p2) < 0; return Cross ( p2 - p1, p3 - p2) < 0;
} }
inline int TRIANGLE2D :: CCW () const inline int TRIANGLE2D :: CCW () const
{ {
return Cross ( p2 - p1, p3 - p2) > 0; return Cross ( p2 - p1, p3 - p2) > 0;
} }
@ -767,19 +768,19 @@ inline int TRIANGLE2D :: CCW () const
inline int PTRIANGLE2D :: Regular() const inline int PTRIANGLE2D :: Regular() const
{ {
return fabs(Cross ( *p2 - *p1, *p3 - *p2)) > EPSGEOM; return fabs(Cross ( *p2 - *p1, *p3 - *p2)) > EPSGEOM;
} }
inline int PTRIANGLE2D :: CW () const inline int PTRIANGLE2D :: CW () const
{ {
return Cross ( *p2 - *p1, *p3 - *p2) < 0; return Cross ( *p2 - *p1, *p3 - *p2) < 0;
} }
inline int PTRIANGLE2D :: CCW () const inline int PTRIANGLE2D :: CCW () const
{ {
return Cross ( *p2 - *p1, *p3 - *p2) > 0; return Cross ( *p2 - *p1, *p3 - *p2) > 0;
} }
@ -788,14 +789,14 @@ inline int PTRIANGLE2D :: CCW () const
#endif #endif
/// ///
class Mat2d class Mat2d
{ {
protected: protected:
/// ///
double coeff[4]; double coeff[4];
public: public:
/// ///
Mat2d() { coeff[0] = coeff[1] = coeff[2] = coeff[3] = 0; } Mat2d() { coeff[0] = coeff[1] = coeff[2] = coeff[3] = 0; }
/// ///
@ -827,27 +828,27 @@ public:
void SolvePositiveDefinite (const Vec2d & rhs, Vec2d & x) const; void SolvePositiveDefinite (const Vec2d & rhs, Vec2d & x) const;
/// add a term \alpha * v * v^T /// add a term \alpha * v * v^T
void AddDiadicProduct (double alpha, Vec2d & v); void AddDiadicProduct (double alpha, Vec2d & v);
}; };
inline void Mat2d :: Mult (const Vec2d & v, Vec2d & prod) const inline void Mat2d :: Mult (const Vec2d & v, Vec2d & prod) const
{ {
prod.X() = coeff[0] * v.X() + coeff[1] * v.Y(); prod.X() = coeff[0] * v.X() + coeff[1] * v.Y();
prod.Y() = coeff[2] * v.X() + coeff[3] * v.Y(); prod.Y() = coeff[2] * v.X() + coeff[3] * v.Y();
} }
inline void Mat2d :: MultTrans (const Vec2d & v, Vec2d & prod) const inline void Mat2d :: MultTrans (const Vec2d & v, Vec2d & prod) const
{ {
prod.X() = coeff[0] * v.X() + coeff[2] * v.Y(); prod.X() = coeff[0] * v.X() + coeff[2] * v.Y();
prod.Y() = coeff[1] * v.X() + coeff[3] * v.Y(); prod.Y() = coeff[1] * v.X() + coeff[3] * v.Y();
} }
inline void Mat2d :: Solve (const Vec2d & rhs, Vec2d & x) const inline void Mat2d :: Solve (const Vec2d & rhs, Vec2d & x) const
{ {
double det = Det(); double det = Det();
if (det == 0) if (det == 0)
@ -857,11 +858,11 @@ inline void Mat2d :: Solve (const Vec2d & rhs, Vec2d & x) const
x.X() = (coeff[3] * rhs.X() - coeff[1] * rhs.Y()) / det; x.X() = (coeff[3] * rhs.X() - coeff[1] * rhs.Y()) / det;
x.Y() = (-coeff[2] * rhs.X() + coeff[0] * rhs.Y()) / det; x.Y() = (-coeff[2] * rhs.X() + coeff[0] * rhs.Y()) / det;
} }
} }
inline void Mat2d :: SolvePositiveDefinite (const Vec2d & rhs, Vec2d & x) const inline void Mat2d :: SolvePositiveDefinite (const Vec2d & rhs, Vec2d & x) const
{ {
double a = max2(coeff[0], 1e-8); double a = max2(coeff[0], 1e-8);
double b = coeff[1] / a; double b = coeff[1] / a;
double c = coeff[2] / a; double c = coeff[2] / a;
@ -869,17 +870,17 @@ inline void Mat2d :: SolvePositiveDefinite (const Vec2d & rhs, Vec2d & x) const
x.X() = (rhs.X() - b * rhs.Y()) / a; x.X() = (rhs.X() - b * rhs.Y()) / a;
x.Y() = rhs.Y() / d - c * x.X(); x.Y() = rhs.Y() / d - c * x.X();
} }
inline void Mat2d :: AddDiadicProduct (double alpha, Vec2d & v) inline void Mat2d :: AddDiadicProduct (double alpha, Vec2d & v)
{ {
coeff[0] += alpha * v.X() * v.X(); coeff[0] += alpha * v.X() * v.X();
coeff[1] += alpha * v.X() * v.Y(); coeff[1] += alpha * v.X() * v.Y();
coeff[2] += alpha * v.Y() * v.X(); coeff[2] += alpha * v.Y() * v.X();
coeff[3] += alpha * v.Y() * v.Y(); coeff[3] += alpha * v.Y() * v.Y();
}
} }
#endif #endif

266
libsrc/gprim/geom3d.hpp
View File

@ -7,69 +7,70 @@
/* Date: 5. Aug. 95 */ /* Date: 5. Aug. 95 */
/* *************************************************************************/ /* *************************************************************************/
namespace netgen
{
extern void MyError (const char * ch);
extern void MyError (const char * ch); class Point3d;
class Vec3d;
class Point3d; inline Vec3d operator- (const Point3d & p1, const Point3d & p2);
class Vec3d; inline Point3d operator- (const Point3d & p1, const Vec3d & v);
inline Point3d operator+ (const Point3d & p1, const Vec3d & v);
inline Vec3d operator- (const Point3d & p1, const Point3d & p2); Point3d & Add (double d, const Vec3d & v);
inline Point3d operator- (const Point3d & p1, const Vec3d & v); Point3d & Add2 (double d, const Vec3d & v,
inline Point3d operator+ (const Point3d & p1, const Vec3d & v);
Point3d & Add (double d, const Vec3d & v);
Point3d & Add2 (double d, const Vec3d & v,
double d2, const Vec3d & v2); double d2, const Vec3d & v2);
inline Point3d Center (const Point3d & p1, const Point3d & p2); inline Point3d Center (const Point3d & p1, const Point3d & p2);
inline Point3d Center (const Point3d & p1, const Point3d & p2, const Point3d & p3); inline Point3d Center (const Point3d & p1, const Point3d & p2, const Point3d & p3);
inline Point3d Center (const Point3d & p1, const Point3d & p2, inline Point3d Center (const Point3d & p1, const Point3d & p2,
const Point3d & p3, const Point3d & p4); const Point3d & p3, const Point3d & p4);
ostream & operator<<(ostream & s, const Point3d & p); ostream & operator<<(ostream & s, const Point3d & p);
inline Vec3d operator- (const Vec3d & p1, const Vec3d & v); inline Vec3d operator- (const Vec3d & p1, const Vec3d & v);
inline Vec3d operator+ (const Vec3d & p1, const Vec3d & v); inline Vec3d operator+ (const Vec3d & p1, const Vec3d & v);
inline Vec3d operator* (double scal, const Vec3d & v); inline Vec3d operator* (double scal, const Vec3d & v);
inline double operator* (const Vec3d & v1, const Vec3d & v2); inline double operator* (const Vec3d & v1, const Vec3d & v2);
inline Vec3d Cross (const Vec3d & v1, const Vec3d & v2); inline Vec3d Cross (const Vec3d & v1, const Vec3d & v2);
inline void Cross (const Vec3d & v1, const Vec3d & v2, Vec3d & prod); inline void Cross (const Vec3d & v1, const Vec3d & v2, Vec3d & prod);
double Angle (const Vec3d & v); double Angle (const Vec3d & v);
double FastAngle (const Vec3d & v); double FastAngle (const Vec3d & v);
double Angle (const Vec3d & v1, const Vec3d & v2); double Angle (const Vec3d & v1, const Vec3d & v2);
double FastAngle (const Vec3d & v1, const Vec3d & v2); double FastAngle (const Vec3d & v1, const Vec3d & v2);
ostream & operator<<(ostream & s, const Vec3d & v); ostream & operator<<(ostream & s, const Vec3d & v);
void Transpose (Vec3d & v1, Vec3d & v2, Vec3d & v3); void Transpose (Vec3d & v1, Vec3d & v2, Vec3d & v3);
int SolveLinearSystem (const Vec3d & col1, int SolveLinearSystem (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec3d & col3, const Vec3d & col3,
const Vec3d & rhs, const Vec3d & rhs,
Vec3d & sol); Vec3d & sol);
int SolveLinearSystemLS (const Vec3d & col1, int SolveLinearSystemLS (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec2d & rhs, const Vec2d & rhs,
Vec3d & sol); Vec3d & sol);
int SolveLinearSystemLS2 (const Vec3d & col1, int SolveLinearSystemLS2 (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec2d & rhs, const Vec2d & rhs,
Vec3d & sol, Vec3d & sol,
double & x, double & y); double & x, double & y);
int PseudoInverse (const Vec3d & col1, int PseudoInverse (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
Vec3d & inv1, Vec3d & inv1,
Vec3d & inv2); Vec3d & inv2);
double Determinant (const Vec3d & col1, double Determinant (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec3d & col3); const Vec3d & col3);
inline double Dist2 (const Point3d & p1, const Point3d & p2); inline double Dist2 (const Point3d & p1, const Point3d & p2);
/// Point in R3 /// Point in R3
class Point3d class Point3d
{ {
protected: protected:
/// ///
double x[3]; double x[3];
public: public:
/// ///
Point3d () { x[0] = x[1] = x[2] = 0; } Point3d () { x[0] = x[1] = x[2] = 0; }
/// ///
@ -176,17 +177,17 @@ public:
{ {
return Point<3> (x[0], x[1], x[2]); return Point<3> (x[0], x[1], x[2]);
} }
}; };
/// ///
class Vec3d class Vec3d
{ {
protected: protected:
/// ///
double x[3]; double x[3];
public: public:
/// ///
inline Vec3d() { x[0] = x[1] = x[2] = 0; } inline Vec3d() { x[0] = x[1] = x[2] = 0; }
/// ///
@ -346,16 +347,16 @@ public:
friend double Determinant (const Vec3d & col1, friend double Determinant (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec3d & col3); const Vec3d & col3);
}; };
class QuadraticFunction3d class QuadraticFunction3d
{ {
double c0, cx, cy, cz; double c0, cx, cy, cz;
double cxx, cyy, czz, cxy, cxz, cyz; double cxx, cyy, czz, cxy, cxz, cyz;
public: public:
QuadraticFunction3d (const Point3d & p, const Vec3d & v); QuadraticFunction3d (const Point3d & p, const Vec3d & v);
double Eval (const Point3d & p) double Eval (const Point3d & p)
{ {
@ -365,64 +366,64 @@ public:
+ p.Y() * (cy + cyy * p.Y() + cyz * p.Z()) + p.Y() * (cy + cyy * p.Y() + cyz * p.Z())
+ p.Z() * (cz + czz * p.Z()); + p.Z() * (cz + czz * p.Z());
} }
}; };
inline Point3d Center (const Point3d & p1, const Point3d & p2) inline Point3d Center (const Point3d & p1, const Point3d & p2)
{ {
return Point3d (0.5 * (p1.x[0] + p2.x[0]), return Point3d (0.5 * (p1.x[0] + p2.x[0]),
0.5 * (p1.x[1] + p2.x[1]), 0.5 * (p1.x[1] + p2.x[1]),
0.5 * (p1.x[2] + p2.x[2])); 0.5 * (p1.x[2] + p2.x[2]));
} }
inline Point3d Center (const Point3d & p1, const Point3d & p2, inline Point3d Center (const Point3d & p1, const Point3d & p2,
const Point3d & p3) const Point3d & p3)
{ {
return Point3d (1.0/3.0 * (p1.x[0] + p2.x[0] + p3.x[0]), return Point3d (1.0/3.0 * (p1.x[0] + p2.x[0] + p3.x[0]),
1.0/3.0 * (p1.x[1] + p2.x[1] + p3.x[1]), 1.0/3.0 * (p1.x[1] + p2.x[1] + p3.x[1]),
1.0/3.0 * (p1.x[2] + p2.x[2] + p3.x[2])); 1.0/3.0 * (p1.x[2] + p2.x[2] + p3.x[2]));
} }
inline Point3d Center (const Point3d & p1, const Point3d & p2, inline Point3d Center (const Point3d & p1, const Point3d & p2,
const Point3d & p3, const Point3d & p4) const Point3d & p3, const Point3d & p4)
{ {
return Point3d (0.25 * (p1.x[0] + p2.x[0] + p3.x[0] + p4.x[0]), return Point3d (0.25 * (p1.x[0] + p2.x[0] + p3.x[0] + p4.x[0]),
0.25 * (p1.x[1] + p2.x[1] + p3.x[1] + p4.x[1]), 0.25 * (p1.x[1] + p2.x[1] + p3.x[1] + p4.x[1]),
0.25 * (p1.x[2] + p2.x[2] + p3.x[2] + p4.x[2])); 0.25 * (p1.x[2] + p2.x[2] + p3.x[2] + p4.x[2]));
} }
inline Vec3d & Vec3d :: operator+= (const Vec3d & v2) inline Vec3d & Vec3d :: operator+= (const Vec3d & v2)
{ {
x[0] += v2.x[0]; x[0] += v2.x[0];
x[1] += v2.x[1]; x[1] += v2.x[1];
x[2] += v2.x[2]; x[2] += v2.x[2];
return *this; return *this;
} }
inline Vec3d & Vec3d :: operator-= (const Vec3d & v2) inline Vec3d & Vec3d :: operator-= (const Vec3d & v2)
{ {
x[0] -= v2.x[0]; x[0] -= v2.x[0];
x[1] -= v2.x[1]; x[1] -= v2.x[1];
x[2] -= v2.x[2]; x[2] -= v2.x[2];
return *this; return *this;
} }
inline Vec3d & Vec3d :: operator*= (double s) inline Vec3d & Vec3d :: operator*= (double s)
{ {
x[0] *= s; x[0] *= s;
x[1] *= s; x[1] *= s;
x[2] *= s; x[2] *= s;
return *this; return *this;
} }
inline Vec3d & Vec3d :: operator/= (double s) inline Vec3d & Vec3d :: operator/= (double s)
{ {
if (s != 0) if (s != 0)
{ {
x[0] /= s; x[0] /= s;
@ -437,24 +438,24 @@ inline Vec3d & Vec3d :: operator/= (double s)
} }
#endif #endif
return *this; return *this;
} }
inline Vec3d & Vec3d::Add (double d, const Vec3d & v) inline Vec3d & Vec3d::Add (double d, const Vec3d & v)
{ {
x[0] += d * v.x[0]; x[0] += d * v.x[0];
x[1] += d * v.x[1]; x[1] += d * v.x[1];
x[2] += d * v.x[2]; x[2] += d * v.x[2];
return *this; return *this;
} }
inline Vec3d & Vec3d::Add2 (double d, const Vec3d & v, inline Vec3d & Vec3d::Add2 (double d, const Vec3d & v,
double d2, const Vec3d & v2) double d2, const Vec3d & v2)
{ {
x[0] += d * v.x[0] + d2 * v2.x[0]; x[0] += d * v.x[0] + d2 * v2.x[0];
x[1] += d * v.x[1] + d2 * v2.x[1]; x[1] += d * v.x[1] + d2 * v2.x[1];
x[2] += d * v.x[2] + d2 * v2.x[2]; x[2] += d * v.x[2] + d2 * v2.x[2];
return *this; return *this;
} }
@ -463,117 +464,117 @@ inline Vec3d & Vec3d::Add2 (double d, const Vec3d & v,
inline Vec3d operator- (const Point3d & p1, const Point3d & p2) inline Vec3d operator- (const Point3d & p1, const Point3d & p2)
{ {
return Vec3d (p1.x[0] - p2.x[0], p1.x[1] - p2.x[1],p1.x[2] - p2.x[2]); return Vec3d (p1.x[0] - p2.x[0], p1.x[1] - p2.x[1],p1.x[2] - p2.x[2]);
} }
inline Point3d operator- (const Point3d & p1, const Vec3d & v) inline Point3d operator- (const Point3d & p1, const Vec3d & v)
{ {
return Point3d (p1.x[0] - v.x[0], p1.x[1] - v.x[1],p1.x[2] - v.x[2]); return Point3d (p1.x[0] - v.x[0], p1.x[1] - v.x[1],p1.x[2] - v.x[2]);
} }
inline Point3d operator+ (const Point3d & p1, const Vec3d & v) inline Point3d operator+ (const Point3d & p1, const Vec3d & v)
{ {
return Point3d (p1.x[0] + v.x[0], p1.x[1] + v.x[1],p1.x[2] + v.x[2]); return Point3d (p1.x[0] + v.x[0], p1.x[1] + v.x[1],p1.x[2] + v.x[2]);
} }
inline Point3d & Point3d::operator+= (const Vec3d & v) inline Point3d & Point3d::operator+= (const Vec3d & v)
{ {
x[0] += v.x[0]; x[0] += v.x[0];
x[1] += v.x[1]; x[1] += v.x[1];
x[2] += v.x[2]; x[2] += v.x[2];
return *this; return *this;
} }
inline Point3d & Point3d::operator-= (const Vec3d & v) inline Point3d & Point3d::operator-= (const Vec3d & v)
{ {
x[0] -= v.x[0]; x[0] -= v.x[0];
x[1] -= v.x[1]; x[1] -= v.x[1];
x[2] -= v.x[2]; x[2] -= v.x[2];
return *this; return *this;
} }
inline Point3d & Point3d::Add (double d, const Vec3d & v) inline Point3d & Point3d::Add (double d, const Vec3d & v)
{ {
x[0] += d * v.x[0]; x[0] += d * v.x[0];
x[1] += d * v.x[1]; x[1] += d * v.x[1];
x[2] += d * v.x[2]; x[2] += d * v.x[2];
return *this; return *this;
} }
inline Point3d & Point3d::Add2 (double d, const Vec3d & v, inline Point3d & Point3d::Add2 (double d, const Vec3d & v,
double d2, const Vec3d & v2) double d2, const Vec3d & v2)
{ {
x[0] += d * v.x[0] + d2 * v2.x[0]; x[0] += d * v.x[0] + d2 * v2.x[0];
x[1] += d * v.x[1] + d2 * v2.x[1]; x[1] += d * v.x[1] + d2 * v2.x[1];
x[2] += d * v.x[2] + d2 * v2.x[2]; x[2] += d * v.x[2] + d2 * v2.x[2];
return *this; return *this;
} }
inline Vec3d operator- (const Vec3d & v1, const Vec3d & v2) inline Vec3d operator- (const Vec3d & v1, const Vec3d & v2)
{ {
return Vec3d (v1.x[0] - v2.x[0], v1.x[1] - v2.x[1],v1.x[2] - v2.x[2]); return Vec3d (v1.x[0] - v2.x[0], v1.x[1] - v2.x[1],v1.x[2] - v2.x[2]);
} }
inline Vec3d operator+ (const Vec3d & v1, const Vec3d & v2) inline Vec3d operator+ (const Vec3d & v1, const Vec3d & v2)
{ {
return Vec3d (v1.x[0] + v2.x[0], v1.x[1] + v2.x[1],v1.x[2] + v2.x[2]); return Vec3d (v1.x[0] + v2.x[0], v1.x[1] + v2.x[1],v1.x[2] + v2.x[2]);
} }
inline Vec3d operator* (double scal, const Vec3d & v) inline Vec3d operator* (double scal, const Vec3d & v)
{ {
return Vec3d (scal * v.x[0], scal * v.x[1], scal * v.x[2]); return Vec3d (scal * v.x[0], scal * v.x[1], scal * v.x[2]);
} }
inline double operator* (const Vec3d & v1, const Vec3d & v2) inline double operator* (const Vec3d & v1, const Vec3d & v2)
{ {
return v1.x[0] * v2.x[0] + v1.x[1] * v2.x[1] + v1.x[2] * v2.x[2]; return v1.x[0] * v2.x[0] + v1.x[1] * v2.x[1] + v1.x[2] * v2.x[2];
} }
inline Vec3d Cross (const Vec3d & v1, const Vec3d & v2) inline Vec3d Cross (const Vec3d & v1, const Vec3d & v2)
{ {
return Vec3d return Vec3d
( v1.x[1] * v2.x[2] - v1.x[2] * v2.x[1], ( v1.x[1] * v2.x[2] - v1.x[2] * v2.x[1],
v1.x[2] * v2.x[0] - v1.x[0] * v2.x[2], v1.x[2] * v2.x[0] - v1.x[0] * v2.x[2],
v1.x[0] * v2.x[1] - v1.x[1] * v2.x[0]); v1.x[0] * v2.x[1] - v1.x[1] * v2.x[0]);
} }
inline void Cross (const Vec3d & v1, const Vec3d & v2, Vec3d & prod) inline void Cross (const Vec3d & v1, const Vec3d & v2, Vec3d & prod)
{ {
prod.x[0] = v1.x[1] * v2.x[2] - v1.x[2] * v2.x[1]; prod.x[0] = v1.x[1] * v2.x[2] - v1.x[2] * v2.x[1];
prod.x[1] = v1.x[2] * v2.x[0] - v1.x[0] * v2.x[2]; prod.x[1] = v1.x[2] * v2.x[0] - v1.x[0] * v2.x[2];
prod.x[2] = v1.x[0] * v2.x[1] - v1.x[1] * v2.x[0]; prod.x[2] = v1.x[0] * v2.x[1] - v1.x[1] * v2.x[0];
} }
inline double Determinant (const Vec3d & col1, inline double Determinant (const Vec3d & col1,
const Vec3d & col2, const Vec3d & col2,
const Vec3d & col3) const Vec3d & col3)
{ {
return return
col1.x[0] * ( col2.x[1] * col3.x[2] - col2.x[2] * col3.x[1]) + col1.x[0] * ( col2.x[1] * col3.x[2] - col2.x[2] * col3.x[1]) +
col1.x[1] * ( col2.x[2] * col3.x[0] - col2.x[0] * col3.x[2]) + col1.x[1] * ( col2.x[2] * col3.x[0] - col2.x[0] * col3.x[2]) +
col1.x[2] * ( col2.x[0] * col3.x[1] - col2.x[1] * col3.x[0]); col1.x[2] * ( col2.x[0] * col3.x[1] - col2.x[1] * col3.x[0]);
} }
/// ///
class Box3d class Box3d
{ {
protected: protected:
/// ///
double minx[3], maxx[3]; double minx[3], maxx[3];
public: public:
/// ///
Box3d () { }; Box3d () { };
/// ///
@ -676,17 +677,17 @@ public:
void WriteData(ofstream& fout) const; void WriteData(ofstream& fout) const;
/// ///
void ReadData(ifstream& fin); void ReadData(ifstream& fin);
}; };
class Box3dSphere : public Box3d class Box3dSphere : public Box3d
{ {
protected: protected:
/// ///
double diam, inner; double diam, inner;
/// ///
Point3d c; Point3d c;
public: public:
/// ///
Box3dSphere () { }; Box3dSphere () { };
/// ///
@ -706,14 +707,14 @@ public:
// private: // private:
/// ///
void CalcDiamCenter (); void CalcDiamCenter ();
}; };
/// ///
class referencetransform class referencetransform
{ {
/// ///
Vec3d ex, ey, ez; Vec3d ex, ey, ez;
/// ///
@ -725,7 +726,7 @@ class referencetransform
/// ///
double h; double h;
public: public:
/// ///
void Set (const Point3d & p1, const Point3d & p2, void Set (const Point3d & p1, const Point3d & p2,
@ -737,8 +738,9 @@ public:
void ToPlain (const Array<Point3d> & p, Array<Point3d> & pp) const; void ToPlain (const Array<Point3d> & p, Array<Point3d> & pp) const;
/// ///
void FromPlain (const Point3d & pp, Point3d & p) const; void FromPlain (const Point3d & pp, Point3d & p) const;
}; };
}
#endif #endif

125
libsrc/gprim/geomfuncs.hpp
View File

@ -8,110 +8,113 @@
/* *************************************************************************/ /* *************************************************************************/
template <int D> namespace netgen
inline double Abs (const Vec<D> & v)
{ {
template <int D>
inline double Abs (const Vec<D> & v)
{
double sum = 0; double sum = 0;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
sum += v(i) * v(i); sum += v(i) * v(i);
return sqrt (sum); return sqrt (sum);
} }
template <int D> template <int D>
inline double Abs2 (const Vec<D> & v) inline double Abs2 (const Vec<D> & v)
{ {
double sum = 0; double sum = 0;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
sum += v(i) * v(i); sum += v(i) * v(i);
return sum; return sum;
} }
template <int D> template <int D>
inline double Dist (const Point<D> & a, const Point<D> & b) inline double Dist (const Point<D> & a, const Point<D> & b)
{ {
return Abs (a-b); return Abs (a-b);
} }
template <int D> template <int D>
inline double Dist2 (const Point<D> & a, const Point<D> & b) inline double Dist2 (const Point<D> & a, const Point<D> & b)
{ {
return Abs2 (a-b); return Abs2 (a-b);
} }
template <int D> template <int D>
inline Point<D> Center (const Point<D> & a, const Point<D> & b) inline Point<D> Center (const Point<D> & a, const Point<D> & b)
{ {
Point<D> res; Point<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = 0.5 * (a(i) + b(i)); res(i) = 0.5 * (a(i) + b(i));
return res; return res;
} }
template <int D> template <int D>
inline Point<D> Center (const Point<D> & a, const Point<D> & b, const Point<D> & c) inline Point<D> Center (const Point<D> & a, const Point<D> & b, const Point<D> & c)
{ {
Point<D> res; Point<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = (1.0/3.0) * (a(i) + b(i) + c(i)); res(i) = (1.0/3.0) * (a(i) + b(i) + c(i));
return res; return res;
} }
template <int D> template <int D>
inline Point<D> Center (const Point<D> & a, const Point<D> & b, const Point<D> & c, const Point<D> & d) inline Point<D> Center (const Point<D> & a, const Point<D> & b, const Point<D> & c, const Point<D> & d)
{ {
Point<D> res; Point<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = (1.0/4.0) * (a(i) + b(i) + c(i) + d(i)); res(i) = (1.0/4.0) * (a(i) + b(i) + c(i) + d(i));
return res; return res;
} }
inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2) inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2)
{ {
return Vec<3> return Vec<3>
( v1(1) * v2(2) - v1(2) * v2(1), ( v1(1) * v2(2) - v1(2) * v2(1),
v1(2) * v2(0) - v1(0) * v2(2), v1(2) * v2(0) - v1(0) * v2(2),
v1(0) * v2(1) - v1(1) * v2(0) ); v1(0) * v2(1) - v1(1) * v2(0) );
} }
inline double Determinant (const Vec<3> & col1, inline double Determinant (const Vec<3> & col1,
const Vec<3> & col2, const Vec<3> & col2,
const Vec<3> & col3) const Vec<3> & col3)
{ {
return return
col1(0) * ( col2(1) * col3(2) - col2(2) * col3(1)) + col1(0) * ( col2(1) * col3(2) - col2(2) * col3(1)) +
col1(1) * ( col2(2) * col3(0) - col2(0) * col3(2)) + col1(1) * ( col2(2) * col3(0) - col2(0) * col3(2)) +
col1(2) * ( col2(0) * col3(1) - col2(1) * col3(0)); col1(2) * ( col2(0) * col3(1) - col2(1) * col3(0));
} }
template <> template <>
inline Vec<2> Vec<2> :: GetNormal () const inline Vec<2> Vec<2> :: GetNormal () const
{ {
return Vec<2> (-x[1], x[0]); return Vec<2> (-x[1], x[0]);
} }
template <> template <>
inline Vec<3> Vec<3> :: GetNormal () const inline Vec<3> Vec<3> :: GetNormal () const
{ {
if (fabs (x[0]) > fabs (x[2])) if (fabs (x[0]) > fabs (x[2]))
return Vec<3> (-x[1], x[0], 0); return Vec<3> (-x[1], x[0], 0);
else else
return Vec<3> (0, x[2], -x[1]); return Vec<3> (0, x[2], -x[1]);
} }
// template <int H, int W> // template <int H, int W>
inline void CalcInverse (const Mat<2,2> & m, Mat<2,2> & inv) inline void CalcInverse (const Mat<2,2> & m, Mat<2,2> & inv)
{ {
double det = m(0,0) * m(1,1) - m(0,1) * m(1,0); double det = m(0,0) * m(1,1) - m(0,1) * m(1,0);
if (det == 0) if (det == 0)
{ {
@ -124,34 +127,36 @@ inline void CalcInverse (const Mat<2,2> & m, Mat<2,2> & inv)
inv(0,1) = -idet * m(0,1); inv(0,1) = -idet * m(0,1);
inv(1,0) = -idet * m(1,0); inv(1,0) = -idet * m(1,0);
inv(1,1) = idet * m(0,0); inv(1,1) = idet * m(0,0);
} }
void CalcInverse (const Mat<3,3> & m, Mat<3,3> & inv); void CalcInverse (const Mat<3,3> & m, Mat<3,3> & inv);
inline void CalcInverse (const Mat<2,3> & m, Mat<3,2> & inv) inline void CalcInverse (const Mat<2,3> & m, Mat<3,2> & inv)
{ {
Mat<2,2> a = m * Trans (m); Mat<2,2> a = m * Trans (m);
Mat<2,2> ainv; Mat<2,2> ainv;
CalcInverse (a, ainv); CalcInverse (a, ainv);
inv = Trans (m) * ainv; inv = Trans (m) * ainv;
} }
void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv); void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv);
inline void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv) inline void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv)
{ {
Mat<2,2> a = Trans (m) * m; Mat<2,2> a = Trans (m) * m;
Mat<2,2> ainv; Mat<2,2> ainv;
CalcInverse (a, ainv); CalcInverse (a, ainv);
inv = ainv * Trans (m); inv = ainv * Trans (m);
}
double Det (const Mat<2,2> & m);
double Det (const Mat<3,3> & m);
// eigenvalues of a symmetric matrix
void EigenValues (const Mat<3,3> & m, Vec<3> & ev);
void EigenValues (const Mat<2,2> & m, Vec<3> & ev);
} }
double Det (const Mat<2,2> & m);
double Det (const Mat<3,3> & m);
// eigenvalues of a symmetric matrix
void EigenValues (const Mat<3,3> & m, Vec<3> & ev);
void EigenValues (const Mat<2,2> & m, Vec<3> & ev);
#endif #endif

73
libsrc/gprim/geomobjects.hpp
View File

@ -8,19 +8,22 @@
/* *************************************************************************/ /* *************************************************************************/
namespace netgen
template <int D> class Vec;
template <int D> class Point;
template <int D>
class Point
{ {
protected:
template <int D> class Vec;
template <int D> class Point;
template <int D>
class Point
{
protected:
double x[D]; double x[D];
public: public:
Point () { ; } Point () { ; }
Point (double ax) { for (int i = 0; i < D; i++) x[i] = ax; } Point (double ax) { for (int i = 0; i < D; i++) x[i] = ax; }
Point (double ax, double ay) { x[0] = ax; x[1] = ay; } Point (double ax, double ay) { x[0] = ax; x[1] = ay; }
@ -52,20 +55,20 @@ public:
const double & operator() (int i) const { return x[i]; } const double & operator() (int i) const { return x[i]; }
operator const double* () const { return x; } operator const double* () const { return x; }
}; };
template <int D> template <int D>
class Vec class Vec
{ {
protected: protected:
double x[D]; double x[D];
public: public:
Vec () { ; } // for (int i = 0; i < D; i++) x[i] = 0; } Vec () { ; } // for (int i = 0; i < D; i++) x[i] = 0; }
Vec (double ax) { for (int i = 0; i < D; i++) x[i] = ax; } Vec (double ax) { for (int i = 0; i < D; i++) x[i] = ax; }
Vec (double ax, double ay) { x[0] = ax; x[1] = ay; } Vec (double ax, double ay) { x[0] = ax; x[1] = ay; }
@ -128,20 +131,20 @@ public:
} }
Vec<D> GetNormal () const; Vec<D> GetNormal () const;
}; };
template <int H, int W=H> template <int H, int W=H>
class Mat class Mat
{ {
protected: protected:
double x[H*W]; double x[H*W];
public: public:
Mat () { ; } Mat () { ; }
Mat (const Mat & b) Mat (const Mat & b)
{ for (int i = 0; i < H*W; i++) x[i] = b.x[i]; } { for (int i = 0; i < H*W; i++) x[i] = b.x[i]; }
@ -185,17 +188,17 @@ public:
CalcInverse (*this, inv); CalcInverse (*this, inv);
sol = inv * rhs; sol = inv * rhs;
} }
}; };
template <int D> template <int D>
class Box class Box
{ {
protected: protected:
Point<D> pmin, pmax; Point<D> pmin, pmax;
public: public:
Box () { ; } Box () { ; }
Box ( const Point<D> & p1) Box ( const Point<D> & p1)
@ -282,22 +285,22 @@ public:
pmax(i) += dist; pmax(i) += dist;
} }
} }
}; };
template <int D> template <int D>
class BoxSphere : public Box<D> class BoxSphere : public Box<D>
{ {
protected: protected:
/// ///
Point<D> c; Point<D> c;
/// ///
double diam; double diam;
/// ///
double inner; double inner;
public: public:
/// ///
BoxSphere () { }; BoxSphere () { };
/// ///
@ -357,11 +360,11 @@ public:
inner = this->pmax(i) - this->pmin(i); inner = this->pmax(i) - this->pmin(i);
} }
}; };
}
#endif #endif

245
libsrc/gprim/geomops.hpp
View File

@ -8,156 +8,159 @@
/* *************************************************************************/ /* *************************************************************************/
/* namespace netgen
{
Point - Vector operations /*
Point - Vector operations
*/ */
template <int D> template <int D>
inline Vec<D> operator+ (const Vec<D> & a, const Vec<D> & b) inline Vec<D> operator+ (const Vec<D> & a, const Vec<D> & b)
{ {
Vec<D> res; Vec<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) + b(i); res(i) = a(i) + b(i);
return res; return res;
} }
template <int D> template <int D>
inline Point<D> operator+ (const Point<D> & a, const Vec<D> & b) inline Point<D> operator+ (const Point<D> & a, const Vec<D> & b)
{ {
Point<D> res; Point<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) + b(i); res(i) = a(i) + b(i);
return res; return res;
} }
template <int D> template <int D>
inline Vec<D> operator- (const Point<D> & a, const Point<D> & b) inline Vec<D> operator- (const Point<D> & a, const Point<D> & b)
{ {
Vec<D> res; Vec<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
} }
template <int D> template <int D>
inline Point<D> operator- (const Point<D> & a, const Vec<D> & b) inline Point<D> operator- (const Point<D> & a, const Vec<D> & b)
{ {
Point<D> res; Point<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
} }
template <int D> template <int D>
inline Vec<D> operator- (const Vec<D> & a, const Vec<D> & b) inline Vec<D> operator- (const Vec<D> & a, const Vec<D> & b)
{ {
Vec<D> res; Vec<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
} }
template <int D> template <int D>
inline Vec<D> operator* (double s, const Vec<D> & b) inline Vec<D> operator* (double s, const Vec<D> & b)
{ {
Vec<D> res; Vec<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = s * b(i); res(i) = s * b(i);
return res; return res;
} }
template <int D> template <int D>
inline double operator* (const Vec<D> & a, const Vec<D> & b) inline double operator* (const Vec<D> & a, const Vec<D> & b)
{ {
double sum = 0; double sum = 0;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
sum += a(i) * b(i); sum += a(i) * b(i);
return sum; return sum;
} }
template <int D> template <int D>
inline Vec<D> operator- (const Vec<D> & b) inline Vec<D> operator- (const Vec<D> & b)
{ {
Vec<D> res; Vec<D> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = -b(i); res(i) = -b(i);
return res; return res;
} }
template <int D> template <int D>
inline Point<D> & operator+= (Point<D> & a, const Vec<D> & b) inline Point<D> & operator+= (Point<D> & a, const Vec<D> & b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) += b(i); a(i) += b(i);
return a; return a;
} }
template <int D> template <int D>
inline Vec<D> & operator+= (Vec<D> & a, const Vec<D> & b) inline Vec<D> & operator+= (Vec<D> & a, const Vec<D> & b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) += b(i); a(i) += b(i);
return a; return a;
} }
template <int D> template <int D>
inline Point<D> & operator-= (Point<D> & a, const Vec<D> & b) inline Point<D> & operator-= (Point<D> & a, const Vec<D> & b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) -= b(i); a(i) -= b(i);
return a; return a;
} }
template <int D> template <int D>
inline Vec<D> & operator-= (Vec<D> & a, const Vec<D> & b) inline Vec<D> & operator-= (Vec<D> & a, const Vec<D> & b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) -= b(i); a(i) -= b(i);
return a; return a;
} }
template <int D> template <int D>
inline Vec<D> & operator*= (Vec<D> & a, double s) inline Vec<D> & operator*= (Vec<D> & a, double s)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) *= s; a(i) *= s;
return a; return a;
} }
template <int D> template <int D>
inline Vec<D> & operator/= (Vec<D> & a, double s) inline Vec<D> & operator/= (Vec<D> & a, double s)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) /= s; a(i) /= s;
return a; return a;
} }
// Matrix - Vector operations // Matrix - Vector operations
/* /*
template <int H, int W> template <int H, int W>
inline Vec<H> operator* (const Mat<H,W> & m, const Vec<W> & v) inline Vec<H> operator* (const Mat<H,W> & m, const Vec<W> & v)
{ {
Vec<H> res; Vec<H> res;
for (int i = 0; i < H; i++) for (int i = 0; i < H; i++)
{ {
@ -166,13 +169,13 @@ inline Vec<H> operator* (const Mat<H,W> & m, const Vec<W> & v)
res(i) += m(i,j) * v(j); res(i) += m(i,j) * v(j);
} }
return res; return res;
} }
*/ */
// thanks to VC60 partial template specialization features !!! // thanks to VC60 partial template specialization features !!!
inline Vec<2> operator* (const Mat<2,2> & m, const Vec<2> & v) inline Vec<2> operator* (const Mat<2,2> & m, const Vec<2> & v)
{ {
Vec<2> res; Vec<2> res;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
{ {
@ -181,10 +184,10 @@ inline Vec<2> operator* (const Mat<2,2> & m, const Vec<2> & v)
res(i) += m(i,j) * v(j); res(i) += m(i,j) * v(j);
} }
return res; return res;
} }
inline Vec<2> operator* (const Mat<2,3> & m, const Vec<3> & v) inline Vec<2> operator* (const Mat<2,3> & m, const Vec<3> & v)
{ {
Vec<2> res; Vec<2> res;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
{ {
@ -193,11 +196,11 @@ inline Vec<2> operator* (const Mat<2,3> & m, const Vec<3> & v)
res(i) += m(i,j) * v(j); res(i) += m(i,j) * v(j);
} }
return res; return res;
} }
inline Vec<3> operator* (const Mat<3,2> & m, const Vec<2> & v) inline Vec<3> operator* (const Mat<3,2> & m, const Vec<2> & v)
{ {
Vec<3> res; Vec<3> res;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
{ {
@ -206,11 +209,11 @@ inline Vec<3> operator* (const Mat<3,2> & m, const Vec<2> & v)
res(i) += m(i,j) * v(j); res(i) += m(i,j) * v(j);
} }
return res; return res;
} }
inline Vec<3> operator* (const Mat<3,3> & m, const Vec<3> & v) inline Vec<3> operator* (const Mat<3,3> & m, const Vec<3> & v)
{ {
Vec<3> res; Vec<3> res;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
{ {
@ -219,7 +222,7 @@ inline Vec<3> operator* (const Mat<3,3> & m, const Vec<3> & v)
res(i) += m(i,j) * v(j); res(i) += m(i,j) * v(j);
} }
return res; return res;
} }
@ -227,10 +230,10 @@ inline Vec<3> operator* (const Mat<3,3> & m, const Vec<3> & v)
/* /*
template <int H1, int W1, int H2, int W2> template <int H1, int W1, int H2, int W2>
inline Mat<H1,W2> operator* (const Mat<H1,W1> & a, const Mat<H2,W2> & b) inline Mat<H1,W2> operator* (const Mat<H1,W1> & a, const Mat<H2,W2> & b)
{ {
Mat<H1,W2> m; Mat<H1,W2> m;
for (int i = 0; i < H1; i++) for (int i = 0; i < H1; i++)
for (int j = 0; j < W2; j++) for (int j = 0; j < W2; j++)
@ -241,11 +244,11 @@ inline Mat<H1,W2> operator* (const Mat<H1,W1> & a, const Mat<H2,W2> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
*/ */
inline Mat<2,2> operator* (const Mat<2,2> & a, const Mat<2,2> & b) inline Mat<2,2> operator* (const Mat<2,2> & a, const Mat<2,2> & b)
{ {
Mat<2,2> m; Mat<2,2> m;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
@ -256,10 +259,10 @@ inline Mat<2,2> operator* (const Mat<2,2> & a, const Mat<2,2> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
inline Mat<2,2> operator* (const Mat<2,3> & a, const Mat<3,2> & b) inline Mat<2,2> operator* (const Mat<2,3> & a, const Mat<3,2> & b)
{ {
Mat<2,2> m; Mat<2,2> m;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
@ -270,11 +273,11 @@ inline Mat<2,2> operator* (const Mat<2,3> & a, const Mat<3,2> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
inline Mat<3,2> operator* (const Mat<3,2> & a, const Mat<2,2> & b) inline Mat<3,2> operator* (const Mat<3,2> & a, const Mat<2,2> & b)
{ {
Mat<3,2> m; Mat<3,2> m;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
@ -285,12 +288,12 @@ inline Mat<3,2> operator* (const Mat<3,2> & a, const Mat<2,2> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
inline Mat<2,3> operator* (const Mat<2,2> & a, const Mat<2,3> & b) inline Mat<2,3> operator* (const Mat<2,2> & a, const Mat<2,3> & b)
{ {
Mat<2,3> m; Mat<2,3> m;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
@ -301,11 +304,11 @@ inline Mat<2,3> operator* (const Mat<2,2> & a, const Mat<2,3> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
inline Mat<3,3> operator* (const Mat<3,3> & a, const Mat<3,3> & b) inline Mat<3,3> operator* (const Mat<3,3> & a, const Mat<3,3> & b)
{ {
Mat<3,3> m; Mat<3,3> m;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
@ -316,7 +319,7 @@ inline Mat<3,3> operator* (const Mat<3,3> & a, const Mat<3,3> & b)
m(i,j) = sum; m(i,j) = sum;
} }
return m; return m;
} }
@ -325,15 +328,15 @@ inline Mat<3,3> operator* (const Mat<3,3> & a, const Mat<3,3> & b)
template <int H, int W> template <int H, int W>
inline Mat<W,H> Trans (const Mat<H,W> & m) inline Mat<W,H> Trans (const Mat<H,W> & m)
{ {
Mat<W,H> res; Mat<W,H> res;
for (int i = 0; i < H; i++) for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++) for (int j = 0; j < W; j++)
res(j,i) = m(i,j); res(j,i) = m(i,j);
return res; return res;
} }
@ -345,36 +348,36 @@ inline Mat<W,H> Trans (const Mat<H,W> & m)
template <int D> template <int D>
inline ostream & operator<< (ostream & ost, const Vec<D> & a) inline ostream & operator<< (ostream & ost, const Vec<D> & a)
{ {
ost << "("; ost << "(";
for (int i = 0; i < D-1; i++) for (int i = 0; i < D-1; i++)
ost << a(i) << ", "; ost << a(i) << ", ";
ost << a(D-1) << ")"; ost << a(D-1) << ")";
return ost; return ost;
} }
template <int D> template <int D>
inline ostream & operator<< (ostream & ost, const Point<D> & a) inline ostream & operator<< (ostream & ost, const Point<D> & a)
{ {
ost << "("; ost << "(";
for (int i = 0; i < D-1; i++) for (int i = 0; i < D-1; i++)
ost << a(i) << ", "; ost << a(i) << ", ";
ost << a(D-1) << ")"; ost << a(D-1) << ")";
return ost; return ost;
} }
template <int D> template <int D>
inline ostream & operator<< (ostream & ost, const Box<D> & b) inline ostream & operator<< (ostream & ost, const Box<D> & b)
{ {
ost << b.PMin() << " - " << b.PMax(); ost << b.PMin() << " - " << b.PMax();
return ost; return ost;
} }
template <int H, int W> template <int H, int W>
inline ostream & operator<< (ostream & ost, const Mat<H,W> & m) inline ostream & operator<< (ostream & ost, const Mat<H,W> & m)
{ {
ost << "("; ost << "(";
for (int i = 0; i < H; i++) for (int i = 0; i < H; i++)
{ {
@ -383,9 +386,9 @@ inline ostream & operator<< (ostream & ost, const Mat<H,W> & m)
ost << endl; ost << endl;
} }
return ost; return ost;
}
} }
#endif #endif

4
libsrc/gprim/geomtest3d.hpp
View File

@ -8,6 +8,9 @@
/* *************************************************************************/ /* *************************************************************************/
namespace netgen
{
extern int extern int
IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line); IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line);
@ -79,5 +82,6 @@ extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2,
extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
const Point3d & l2p1, const Point3d & l2p2); const Point3d & l2p1, const Point3d & l2p2);
}
#endif #endif

17
libsrc/gprim/gprim.hpp
View File

@ -8,19 +8,26 @@
/* *************************************************************************/ /* *************************************************************************/
namespace netgen #include <myadt.hpp>
{ #include <linalg.hpp>
#include "geomobjects.hpp" #include "geomobjects.hpp"
#include "geomops.hpp" #include "geomops.hpp"
#include "geomfuncs.hpp" #include "geomfuncs.hpp"
#include "geom2d.hpp" #include "geom2d.hpp"
#include "geom3d.hpp" #include "geom3d.hpp"
#include "geomtest3d.hpp" #include "geomtest3d.hpp"
// #include "rot3d.hpp"
#include "transform3d.hpp" #include "transform3d.hpp"
// #include "reftrans.hpp"
#include "adtree.hpp" #include "adtree.hpp"
}
#include "spline.hpp"
#include "splinegeometry.hpp"
#endif #endif

5
libsrc/gprim/transform3d.hpp
View File

@ -11,6 +11,9 @@
Affine - Linear mapping in 3D space Affine - Linear mapping in 3D space
*/ */
namespace netgen
{
class Transformation3d; class Transformation3d;
ostream & operator<< (ostream & ost, Transformation3d & trans); ostream & operator<< (ostream & ost, Transformation3d & trans);
@ -185,6 +188,6 @@ template <int D>
ostream & operator<< (ostream & ost, Transformation<D> & trans); ostream & operator<< (ostream & ost, Transformation<D> & trans);
}
#endif #endif

11
libsrc/stlgeom/stlgeom.cpp
View File

@ -1,14 +1,7 @@
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <gprim.hpp>
#include <meshing.hpp> #include <meshing.hpp>
#include "stlgeom.hpp" #include "stlgeom.hpp"
namespace netgen namespace netgen
{ {
@ -39,12 +32,14 @@ void STLMeshing (STLGeometry & geom,
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
STLGeometry :: STLGeometry() STLGeometry :: STLGeometry()
/*
: edges(), edgesperpoint(), : edges(), edgesperpoint(),
normals(), externaledges(), normals(), externaledges(),
atlas(), chartmark(), atlas(), chartmark(),
lines(), outerchartspertrig(), vicinity(), markedtrigs(), markedsegs(), lines(), outerchartspertrig(), vicinity(), markedtrigs(), markedsegs(),
lineendpoints(), spiralpoints(), selectedmultiedge() lineendpoints(), spiralpoints(), selectedmultiedge()
*/
{ {
edgedata = new STLEdgeDataList(*this); edgedata = new STLEdgeDataList(*this);
externaledges.SetSize(0); externaledges.SetSize(0);

2
libsrc/stlgeom/stlgeom.hpp
View File

@ -21,11 +21,9 @@
*/ */
#include <gprim.hpp>
#include <meshing.hpp> #include <meshing.hpp>
namespace netgen namespace netgen
{ {
extern int IsInArray(int n, const Array<int>& ia); extern int IsInArray(int n, const Array<int>& ia);

3
libsrc/stlgeom/stltool.cpp
View File

@ -856,7 +856,8 @@ void STLBoundarySeg :: Swap ()
STLBoundary :: STLBoundary (STLGeometry * ageometry) STLBoundary :: STLBoundary (STLGeometry * ageometry)
: boundary(), geometry(ageometry) : // boundary(),
geometry(ageometry)
{ {
; ;
} }