Skip to content
Snippets Groups Projects
Commit 4f157b3c authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent fd733494
No related branches found
No related tags found
No related merge requests found
Showing
with 0 additions and 12459 deletions
#
#
appl = NETGEN
.default all:
#
#
all:
@ (cd linalg; $(MAKE) -f Makefile)
@ (cd general; $(MAKE) -f Makefile)
@ (cd gprim; $(MAKE) -f Makefile)
@ (cd csg; $(MAKE) -f Makefile)
@ (cd geom2d; $(MAKE) -f Makefile)
@ (cd stlgeom; $(MAKE) -f Makefile)
@ (cd occ; $(MAKE) -f Makefile)
@ (cd meshing; $(MAKE) -f Makefile)
@ (cd opti; $(MAKE) -f Makefile)
@ (cd visualization; $(MAKE) -f Makefile)
@ (cd interface; $(MAKE) -f Makefile)
#
# @ (cd step; $(MAKE) -f Makefile)
# @ (cd stepgeom; $(MAKE) -f Makefile)
# @ (cd graphics; $(MAKE) -f Makefile)
tar:
tar cvf ../../libsrc.tar Makefile
tar rf ../../libsrc.tar linalg/Makefile linalg/*.hh linalg/*.cc
tar rf ../../libsrc.tar general/Makefile general/*.hh general/*.cc
tar rf ../../libsrc.tar gprim/Makefile gprim/*.hh gprim/*.cc
tar rf ../../libsrc.tar csg/Makefile csg/*.hh csg/*.cc
tar rf ../../libsrc.tar stlgeom/Makefile stlgeom/*.hh stlgeom/*.cc
tar rf ../../libsrc.tar occ/Makefile occ/*.h* occ/*.c*
tar rf ../../libsrc.tar meshing/Makefile meshing/*.hh meshing/*.cc meshing/*.h
tar rf ../../libsrc.tar opti/Makefile opti/*.hh opti/*.cc
tar rf ../../libsrc.tar step/Makefile step/*.h step/*.cc
tar rf ../../libsrc.tar stepgeom/Makefile stepgeom/*.hh stepgeom/*.cc
tar tf ../../libsrc.tar include/*.h include/*.hh
gzip -9 ../../libsrc.tar
#
# Makefile for geometric library
#
src = csgparser.cpp algprim.cpp curve2d.cpp brick.cpp \
solid.cpp spline3d.cpp surface.cpp bspline2d.cpp \
explicitcurve2d.cpp gencyl.cpp csgeom.cpp polyhedra.cpp extrusion.cpp revolution.cpp \
manifold.cpp curve2d.cpp triapprox.cpp identify.cpp \
singularref.cpp \
edgeflw.cpp specpoin.cpp meshsurf.cpp genmesh.cpp
#
# lex.yy.cpp geometry.cpp
#
lib = csg
libpath = libsrc/csg
#
#
include ../makefile.inc
#
# geometry.cpp : geometry.yy
# bison -d -o geometry.c geometry.yy
# mv -f geometry.c geometry.cpp
#
# lex.yy.cpp : geometry.yy geometry.ll
# flex -+ -d -I geometry.ll
# mv lex.yy.cc lex.yy.cpp
This diff is collapsed.
#ifndef FILE_ALGPRIM
#define FILE_ALGPRIM
/**************************************************************************/
/* File: algprim.hh */
/* Author: Joachim Schoeberl */
/* Date: 1. Dez. 95 */
/**************************************************************************/
/*
Quadric Surfaces (Plane, Sphere, Cylinder)
*/
/**
A quadric surface.
surface defined by
cxx x^2 + cyy y^2 + czz z^2 + cxy x y + cxz x z + cyz y z +
cx x + cy y + cz z + c1 = 0.
**/
class QuadraticSurface : public OneSurfacePrimitive
{
protected:
double cxx, cyy, czz, cxy, cxz, cyz, cx, cy, cz, c1;
public:
virtual double CalcFunctionValue (const Point<3> & point) const;
virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const;
virtual void CalcHesse (const Point<3> & point, Mat<3> & hesse) const;
/*
virtual int RootInBox (const Box<3> & box)
const { return 0; }
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box)
const { return DOES_INTERSECT; }
*/
virtual double HesseNorm () const { return cxx + cyy + czz; }
virtual Point<3> GetSurfacePoint () const;
virtual void Print (ostream & ist) const;
virtual void Read (istream & ist);
void PrintCoeff (ostream & ost) const;
};
/// A Plane (i.e., the plane and everything behind it).
class Plane : public QuadraticSurface
{
/// a point in the plane
Point<3> p;
/// outward normal vector
Vec<3> n;
public:
///
Plane (const Point<3> & ap, Vec<3> an);
virtual void GetPrimitiveData (char *& classname,
ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
static Primitive * CreateDefault ();
virtual Primitive * Copy () const;
virtual void Transform (Transformation<3> & trans);
virtual int IsIdentic (const Surface & s2, int & inv, double eps) const;
///
virtual void DefineTangentialPlane (const Point<3> & ap1,
const Point<3> & ap2);
///
virtual void ToPlane (const Point<3> & p3d,
Point<2> & pplane, double h,
int & zone) const;
///
virtual void FromPlane (const Point<2> & pplane,
Point<3> & p3d,
double h) const;
///
virtual void Project (Point<3> & p) const;
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
inline virtual double CalcFunctionValue (const Point<3> & p3d) const
{return cx * p3d(0) + cy * p3d(1) + cz * p3d(2) + c1;}
///
virtual void CalcGradient (const Point<3> & point,
Vec<3> & grad) const;
///
virtual void CalcHesse (const Point<3> & point,
Mat<3> & hesse) const;
///
virtual double HesseNorm () const;
///
virtual Point<3> GetSurfacePoint () const;
///
virtual void GetTriangleApproximation
(TriangleApproximation & tas,
const Box<3> & boundingbox, double facets) const;
};
// typedef Plane Plane;
///
class Sphere : public QuadraticSurface
{
///
Point<3> c;
///
double r;
public:
///
Sphere (const Point<3> & ac, double ar);
virtual void GetPrimitiveData (char *& classname,
ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
static Primitive * CreateDefault ();
virtual Primitive * Copy () const;
virtual void Transform (Transformation<3> & trans);
virtual int IsIdentic (const Surface & s2, int & inv, double eps) const;
///
virtual void DefineTangentialPlane (const Point<3> & ap1,
const Point<3> & ap2);
///
virtual void ToPlane (const Point<3> & p3d,
Point<2> & pplane, double h,
int & zone) const;
///
virtual void FromPlane (const Point<2> & pplane,
Point<3> & p, double h) const;
///
virtual void Project (Point<3> & p) const;
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
virtual double HesseNorm () const;
///
virtual Point<3> GetSurfacePoint () const;
///
const Point<3> & Center () const { return c; }
///
double Radius () const { return r; }
///
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const;
};
///
class Cylinder : public QuadraticSurface
{
///
Point<3> a, b;
///
double r;
///
Vec<3> vab;
public:
Cylinder (const Point<3> & aa, const Point<3> & ab, double ar);
virtual void GetPrimitiveData (char *& classname, ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
static Primitive * CreateDefault ();
virtual Primitive * Copy () const;
virtual void Transform (Transformation<3> & trans);
///
virtual int IsIdentic (const Surface & s2, int & inv, double eps) const;
///
virtual void DefineTangentialPlane (const Point<3> & ap1,
const Point<3> & ap2);
///
virtual void ToPlane (const Point<3> & p,
Point<2> & pplane,
double h,
int & zone) const;
///
virtual void FromPlane (const Point<2> & pplane,
Point<3> & p,
double h) const;
///
virtual void Project (Point<3> & p) const;
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
virtual double HesseNorm () const;
///
virtual Point<3> GetSurfacePoint () const;
///
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const;
};
///
class EllipticCylinder : public QuadraticSurface
{
private:
///
Point<3> a;
///
Vec<3> vl, vs;
///
Vec<3> vab, t0vec, t1vec;
///
double vabl, t0, t1;
public:
///
EllipticCylinder (const Point<3> & aa,
const Vec<3> & avl, const Vec<3> & avs);
/*
static Primitive * CreateDefault ();
virtual void GetPrimitiveData (char *& classname, ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
*/
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
virtual double HesseNorm () const;
///
virtual Point<3> GetSurfacePoint () const;
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const;
private:
void CalcData();
};
///
class Ellipsoid : public QuadraticSurface
{
private:
///
Point<3> a;
///
Vec<3> v1, v2, v3;
///
double rmin;
public:
///
Ellipsoid (const Point<3> & aa,
const Vec<3> & av1,
const Vec<3> & av2,
const Vec<3> & av3);
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
virtual double HesseNorm () const;
///
virtual Point<3> GetSurfacePoint () const;
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const;
private:
void CalcData();
};
///
class Cone : public QuadraticSurface
{
///
Point<3> a, b;
///
double ra, rb, minr;
///
Vec<3> vab, t0vec, t1vec;
///
double vabl, t0, t1;
public:
///
Cone (const Point<3> & aa, const Point<3> & ab, double ara, double arb);
///
static Primitive * CreateDefault ();
virtual void GetPrimitiveData (char *& classname, ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
///
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
///
virtual double HesseNorm () const;
virtual double LocH (const Point<3> & p, double x,
double c, double hmax) const;
///
virtual Point<3> GetSurfacePoint () const;
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const;
private:
void CalcData();
};
#endif
#include <mystdlib.h>
#include <linalg.hpp>
#include <csg.hpp>
namespace netgen
{
Parallelogram3d :: Parallelogram3d (Point<3> ap1, Point<3> ap2, Point<3> ap3)
{
p1 = ap1;
p2 = ap2;
p3 = ap3;
CalcData();
}
Parallelogram3d ::~Parallelogram3d ()
{
;
}
void Parallelogram3d :: SetPoints (Point<3> ap1,
Point<3> ap2,
Point<3> ap3)
{
p1 = ap1;
p2 = ap2;
p3 = ap3;
CalcData();
}
void Parallelogram3d :: CalcData()
{
v12 = p2 - p1;
v13 = p3 - p1;
p4 = p2 + v13;
n = Cross (v12, v13);
n.Normalize();
}
int Parallelogram3d ::
IsIdentic (const Surface & s2, int & inv, double eps) const
{
int id =
(fabs (s2.CalcFunctionValue (p1)) <= eps) &&
(fabs (s2.CalcFunctionValue (p2)) <= eps) &&
(fabs (s2.CalcFunctionValue (p3)) <= eps);
if (id)
{
Vec<3> n2;
n2 = s2.GetNormalVector(p1);
inv = (n * n2) < 0;
}
return id;
}
double Parallelogram3d :: CalcFunctionValue (const Point<3> & point) const
{
return n * (point - p1);
}
void Parallelogram3d :: CalcGradient (const Point<3> & point,
Vec<3> & grad) const
{
grad = n;
}
void Parallelogram3d :: CalcHesse (const Point<3> & point, Mat<3> & hesse) const
{
hesse = 0;
}
double Parallelogram3d :: HesseNorm () const
{
return 0;
}
Point<3> Parallelogram3d :: GetSurfacePoint () const
{
return p1;
}
void Parallelogram3d :: Print (ostream & str) const
{
str << "Parallelogram3d " << p1 << " - " << p2 << " - " << p3 << endl;
}
void Parallelogram3d ::
GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & bbox,
double facets) const
{
tas.AddPoint (p1);
tas.AddPoint (p2);
tas.AddPoint (p3);
tas.AddPoint (p4);
tas.AddTriangle (TATriangle (0, 0, 1, 2));
tas.AddTriangle (TATriangle (0, 2, 1, 3));
}
Brick :: Brick (Point<3> ap1, Point<3> ap2,
Point<3> ap3, Point<3> ap4)
{
faces.SetSize (6);
surfaceids.SetSize (6);
surfaceactive.SetSize(6);
p1 = ap1; p2 = ap2;
p3 = ap3; p4 = ap4;
for (int i = 0; i < 6; i++)
{
faces[i] = new Plane (Point<3>(0,0,0), Vec<3> (0,0,1));
surfaceactive[i] = 1;
}
CalcData();
}
Brick :: ~Brick ()
{
for (int i = 0; i < 6; i++)
delete faces[i];
}
Primitive * Brick :: CreateDefault ()
{
return new Brick (Point<3> (0,0,0),
Point<3> (1,0,0),
Point<3> (0,1,0),
Point<3> (0,0,1));
}
INSOLID_TYPE Brick :: BoxInSolid (const BoxSphere<3> & box) const
{
/*
int i;
double maxval;
for (i = 1; i <= 6; i++)
{
double val = faces.Get(i)->CalcFunctionValue (box.Center());
if (i == 1 || val > maxval)
maxval = val;
}
if (maxval > box.Diam()) return IS_OUTSIDE;
if (maxval < -box.Diam()) return IS_INSIDE;
return DOES_INTERSECT;
*/
bool inside = 1;
bool outside = 0;
for (int i = 0; i < 6; i++)
{
bool outsidei = 1;
for (int j = 0; j < 8; j++)
{
Point<3> p = box.GetPointNr (j);
double val = faces[i]->CalcFunctionValue (p);
if (val > 0) inside = 0;
if (val < 0) outsidei = 0;
}
if (outsidei) outside = 1;
}
if (outside) return IS_OUTSIDE;
if (inside) return IS_INSIDE;
return DOES_INTERSECT;
}
INSOLID_TYPE Brick :: PointInSolid (const Point<3> & p,
double eps) const
{
double maxval = faces[0] -> CalcFunctionValue (p);
for (int i = 1; i < 6; i++)
{
double val = faces[i] -> CalcFunctionValue (p);
if (val > maxval) maxval = val;
}
if (maxval > eps) return IS_OUTSIDE;
if (maxval < -eps) return IS_INSIDE;
return DOES_INTERSECT;
}
INSOLID_TYPE Brick :: VecInSolid (const Point<3> & p,
const Vec<3> & v,
double eps) const
{
INSOLID_TYPE is = IS_INSIDE;
Vec<3> grad;
double scal;
for (int i = 0; i < faces.Size(); i++)
{
if (faces[i] -> PointOnSurface (p, eps))
{
GetSurface(i).CalcGradient (p, grad);
scal = v * grad;
if (scal >= eps)
is = IS_OUTSIDE;
if (scal >= -eps && is == IS_INSIDE)
is = DOES_INTERSECT;
}
}
return is;
/*
Point<3> p2 = p + 1e-2 * v;
return PointInSolid (p2, eps);
*/
}
void Brick ::
GetPrimitiveData (char *& classname, ARRAY<double> & coeffs) const
{
classname = "brick";
coeffs.SetSize(12);
coeffs.Elem(1) = p1(0);
coeffs.Elem(2) = p1(1);
coeffs.Elem(3) = p1(2);
coeffs.Elem(4) = p2(0);
coeffs.Elem(5) = p2(1);
coeffs.Elem(6) = p2(2);
coeffs.Elem(7) = p3(0);
coeffs.Elem(8) = p3(1);
coeffs.Elem(9) = p3(2);
coeffs.Elem(10) = p4(0);
coeffs.Elem(11) = p4(1);
coeffs.Elem(12) = p4(2);
}
void Brick :: SetPrimitiveData (ARRAY<double> & coeffs)
{
p1(0) = coeffs.Elem(1);
p1(1) = coeffs.Elem(2);
p1(2) = coeffs.Elem(3);
p2(0) = coeffs.Elem(4);
p2(1) = coeffs.Elem(5);
p2(2) = coeffs.Elem(6);
p3(0) = coeffs.Elem(7);
p3(1) = coeffs.Elem(8);
p3(2) = coeffs.Elem(9);
p4(0) = coeffs.Elem(10);
p4(1) = coeffs.Elem(11);
p4(2) = coeffs.Elem(12);
CalcData();
}
void Brick :: CalcData()
{
v12 = p2 - p1;
v13 = p3 - p1;
v14 = p4 - p1;
Point<3> pi[8];
int i1, i2, i3;
int i, j;
i = 0;
for (i3 = 0; i3 <= 1; i3++)
for (i2 = 0; i2 <= 1; i2++)
for (i1 = 0; i1 <= 1; i1++)
{
pi[i] = p1 + i1 * v12 + i2 * v13 + i3 * v14;
i++;
}
static int lface[6][4] =
{ { 1, 3, 2, 4 },
{ 5, 6, 7, 8 },
{ 1, 2, 5, 6 },
{ 3, 7, 4, 8 },
{ 1, 5, 3, 7 },
{ 2, 4, 6, 8 } };
ARRAY<double> data(6);
for (i = 0; i < 6; i++)
{
const Point<3> p1 = pi[lface[i][0]-1];
const Point<3> p2 = pi[lface[i][1]-1];
const Point<3> p3 = pi[lface[i][2]-1];
Vec<3> n = Cross ((p2-p1), (p3-p1));
n.Normalize();
for (j = 0; j < 3; j++)
{
data[j] = p1(j);
data[j+3] = n(j);
}
faces[i] -> SetPrimitiveData (data);
/*
{
faces.Elem(i+1) -> SetPoints
(pi[lface[i][0]-1],
pi[lface[i][1]-1],
pi[lface[i][2]-1]);
}
*/
}
}
void Brick :: Reduce (const BoxSphere<3> & box)
{
double val;
Point<3> p;
for (int i = 0; i < 6; i++)
{
bool hasout = 0;
bool hasin = 0;
for (int j = 0; j < 8; j++)
{
p = box.GetPointNr (j);
val = faces[i]->CalcFunctionValue (p);
if (val > 0) hasout = 1;
else if (val < 0) hasin = 1;
}
surfaceactive[i] = hasout && hasin;
}
}
void Brick :: UnReduce ()
{
for (int i = 0; i < 6; i++)
surfaceactive[i] = 1;
}
OrthoBrick :: OrthoBrick (const Point<3> & ap1, const Point<3> & ap2)
: Brick (ap1,
Point<3> (ap2(0), ap1(1), ap1(2)),
Point<3> (ap1(0), ap2(1), ap1(2)),
Point<3> (ap1(0), ap1(1), ap2(2)))
{
pmin = ap1;
pmax = ap2;
}
INSOLID_TYPE OrthoBrick :: BoxInSolid (const BoxSphere<3> & box) const
{
if (pmin(0) > box.PMax()(0) ||
pmin(1) > box.PMax()(1) ||
pmin(2) > box.PMax()(2) ||
pmax(0) < box.PMin()(0) ||
pmax(1) < box.PMin()(1) ||
pmax(2) < box.PMin()(2))
return IS_OUTSIDE;
if (pmin(0) < box.PMin()(0) &&
pmin(1) < box.PMin()(1) &&
pmin(2) < box.PMin()(2) &&
pmax(0) > box.PMax()(0) &&
pmax(1) > box.PMax()(1) &&
pmax(2) > box.PMax()(2))
return IS_INSIDE;
return DOES_INTERSECT;
}
void OrthoBrick :: Reduce (const BoxSphere<3> & box)
{
surfaceactive.Elem(1) =
(box.PMin()(2) < pmin(2)) && (pmin(2) < box.PMax()(2));
surfaceactive.Elem(2) =
(box.PMin()(2) < pmax(2)) && (pmax(2) < box.PMax()(2));
surfaceactive.Elem(3) =
(box.PMin()(1) < pmin(1)) && (pmin(1) < box.PMax()(1));
surfaceactive.Elem(4) =
(box.PMin()(1) < pmax(1)) && (pmax(1) < box.PMax()(1));
surfaceactive.Elem(5) =
(box.PMin()(0) < pmin(0)) && (pmin(0) < box.PMax()(0));
surfaceactive.Elem(6) =
(box.PMin()(0) < pmax(0)) && (pmax(0) < box.PMax()(0));
}
}
#ifndef FILE_BRICK
#define FILE_BRICK
/**************************************************************************/
/* File: brick.hh */
/* Author: Joachim Schoeberl */
/* Date: 11. Mar. 98 */
/**************************************************************************/
/*
brick geometry, has several surfaces
*/
class Parallelogram3d : public Surface
{
Point<3> p1, p2, p3, p4;
Vec<3> v12, v13;
Vec<3> n;
public:
Parallelogram3d (Point<3> ap1, Point<3> ap2, Point<3> ap3);
virtual ~Parallelogram3d ();
void SetPoints (Point<3> ap1, Point<3> ap2, Point<3> ap3);
virtual int IsIdentic (const Surface & s2, int & inv, double eps) const;
virtual double CalcFunctionValue (const Point<3> & point) const;
virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const;
virtual void CalcHesse (const Point<3> & point, Mat<3> & hesse) const;
virtual double HesseNorm () const;
virtual Point<3> GetSurfacePoint () const;
virtual void Print (ostream & str) const;
virtual void GetTriangleApproximation (TriangleApproximation & tas,
const Box<3> & boundingbox,
double facets) const;
protected:
void CalcData();
};
class Brick : public Primitive
{
Point<3> p1, p2, p3, p4;
Vec<3> v12, v13, v14;
ARRAY<OneSurfacePrimitive*> faces;
public:
Brick (Point<3> ap1, Point<3> ap2, Point<3> ap3, Point<3> ap4);
virtual ~Brick ();
static Primitive * CreateDefault ();
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
virtual INSOLID_TYPE PointInSolid (const Point<3> & p,
double eps) const;
virtual INSOLID_TYPE VecInSolid (const Point<3> & p,
const Vec<3> & v,
double eps) const;
virtual int GetNSurfaces() const
{ return 6; }
virtual Surface & GetSurface (int i)
{ return *faces[i]; }
virtual const Surface & GetSurface (int i) const
{ return *faces[i]; }
virtual void GetPrimitiveData (char *& classname, ARRAY<double> & coeffs) const;
virtual void SetPrimitiveData (ARRAY<double> & coeffs);
virtual void Reduce (const BoxSphere<3> & box);
virtual void UnReduce ();
protected:
void CalcData();
};
class OrthoBrick : public Brick
{
protected:
Point<3> pmin, pmax;
public:
OrthoBrick (const Point<3> & ap1, const Point<3> & ap2);
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
virtual void Reduce (const BoxSphere<3> & box);
};
#endif
#include <mystdlib.h>
#include <csg.hpp>
namespace netgen
{
BSplineCurve2d :: BSplineCurve2d ()
{
redlevel = 0;
}
void BSplineCurve2d :: AddPoint (const Point<2> & apoint)
{
points.Append (apoint);
intervallused.Append (0);
}
bool BSplineCurve2d :: Inside (const Point<2> & p, double & dist) const
{
Point<2> hp = p;
double t = ProjectParam (p);
hp = Eval(t);
Vec<2> v = EvalPrime (t);
Vec<2> n (v(0), -v(1));
cout << "p = " << p << ", hp = " << hp << endl;
dist = Dist (p, hp);
double scal = (hp-p) * n;
cout << "scal = " << scal << endl;
return scal >= 0;
}
double BSplineCurve2d :: ProjectParam (const Point<2> & p) const
{
double t, dt, mindist, mint = 0.0;
int n1;
mindist = 1e10;
dt = 0.2;
for (n1 = 1; n1 <= points.Size(); n1++)
if (intervallused.Get(n1) == 0)
for (t = n1; t <= n1+1; t += dt)
if (Dist (Eval(t), p) < mindist)
{
mint = t;
mindist = Dist (Eval(t), p);
}
if (mindist > 1e9)
{
for (t = 0; t <= points.Size(); t += dt)
if (Dist (Eval(t), p) < mindist)
{
mint = t;
mindist = Dist (Eval(t), p);
}
}
while (Dist (Eval (mint-dt), p) < mindist)
{
mindist = Dist (Eval (mint-dt), p);
mint -= dt;
}
while (Dist (Eval (mint+dt), p) < mindist)
{
mindist = Dist (Eval (mint+dt), p);
mint += dt;
}
return NumericalProjectParam (p, mint-dt, mint+dt);
}
// t \in (n1, n2)
Point<2> BSplineCurve2d :: Eval (double t) const
{
int n, n1, n2, n3, n4;
double loct, b1, b2, b3, b4;
Point<2> hp;
static int cnt = 0;
cnt++;
if (cnt % 100000 == 0) (*mycout) << "cnt = " << cnt << endl;
n = int(t);
loct = t - n;
b1 = 0.25 * (1 - loct) * (1 - loct);
b4 = 0.25 * loct * loct;
b2 = 0.5 - b4;
b3 = 0.5 - b1;
n1 = (n + 10 * points.Size() -1) % points.Size() + 1;
n2 = n1+1;
if (n2 > points.Size()) n2 = 1;
n3 = n2+1;
if (n3 > points.Size()) n3 = 1;
n4 = n3+1;
if (n4 > points.Size()) n4 = 1;
// (*mycout) << "t = " << t << " n = " << n << " loct = " << loct
// << " n1 = " << n1 << endl;
hp(0) = b1 * points.Get(n1)(0) + b2 * points.Get(n2)(0) +
b3 * points.Get(n3)(0) + b4 * points.Get(n4)(0);
hp(1) = b1 * points.Get(n1)(1) + b2 * points.Get(n2)(1) +
b3 * points.Get(n3)(1) + b4 * points.Get(n4)(1);
return hp;
}
Vec<2> BSplineCurve2d :: EvalPrime (double t) const
{
int n, n1, n2, n3, n4;
double loct, db1, db2, db3, db4;
Vec<2> hv;
n = int(t);
loct = t - n;
db1 = 0.5 * (loct - 1);
db4 = 0.5 * loct;
db2 = -db4;
db3 = -db1;
n1 = (n + 10 * points.Size() -1) % points.Size() + 1;
n2 = n1+1;
if (n2 > points.Size()) n2 = 1;
n3 = n2+1;
if (n3 > points.Size()) n3 = 1;
n4 = n3+1;
if (n4 > points.Size()) n4 = 1;
hv(0) = db1 * points.Get(n1)(0) + db2 * points.Get(n2)(0) +
db3 * points.Get(n3)(0) + db4 * points.Get(n4)(0);
hv(1) = db1 * points.Get(n1)(1) + db2 * points.Get(n2)(1) +
db3 * points.Get(n3)(1) + db4 * points.Get(n4)(1);
return hv;
}
Vec<2> BSplineCurve2d :: EvalPrimePrime (double t) const
{
int n, n1, n2, n3, n4;
double ddb1, ddb2, ddb3, ddb4;
Vec<2> hv;
n = int(t);
// double loct = t - n;
ddb1 = 0.5;
ddb4 = 0.5;
ddb2 = -0.5;
ddb3 = -0.5;
n1 = (n + 10 * points.Size() -1) % points.Size() + 1;
n2 = n1+1;
if (n2 > points.Size()) n2 = 1;
n3 = n2+1;
if (n3 > points.Size()) n3 = 1;
n4 = n3+1;
if (n4 > points.Size()) n4 = 1;
hv(0) = ddb1 * points.Get(n1)(0) + ddb2 * points.Get(n2)(0) +
ddb3 * points.Get(n3)(0) + ddb4 * points.Get(n4)(0);
hv(1) = ddb1 * points.Get(n1)(1) + ddb2 * points.Get(n2)(1) +
ddb3 * points.Get(n3)(1) + ddb4 * points.Get(n4)(1);
return hv;
}
int BSplineCurve2d :: SectionUsed (double t) const
{
int n1 = int(t);
n1 = (n1 + 10 * points.Size() - 1) % points.Size() + 1;
return (intervallused.Get(n1) == 0);
}
void BSplineCurve2d :: Reduce (const Point<2> & p, double rad)
{
int n1, n;
int j;
double minx, miny, maxx, maxy;
// (*testout) << "Reduce: " << p << "," << rad << endl;
redlevel++;
for (n1 = 1; n1 <= points.Size(); n1++)
{
if (intervallused.Get(n1) != 0) continue;
minx = maxx = points.Get(n1)(0);
miny = maxy = points.Get(n1)(1);
n = n1;
for (j = 1; j <= 3; j++)
{
n++;
if (n > points.Size()) n = 1;
if (points.Get(n)(0) < minx) minx = points.Get(n)(0);
if (points.Get(n)(1) < miny) miny = points.Get(n)(1);
if (points.Get(n)(0) > maxx) maxx = points.Get(n)(0);
if (points.Get(n)(1) > maxy) maxy = points.Get(n)(1);
}
if (minx > p(0) + rad || maxx < p(0) - rad ||
miny > p(1) + rad || maxy < p(1) - rad)
{
intervallused.Elem(n1) = redlevel;
// (*testout) << 0;
}
else
{
// (*testout) << 1;
intervallused.Elem(n1) = 0;
}
}
// (*testout) << endl;
}
void BSplineCurve2d :: UnReduce ()
{
int i;
for (i = 1; i <= intervallused.Size(); i++)
if (intervallused.Get(i) == redlevel)
intervallused.Set (i, 0);
redlevel--;
}
void BSplineCurve2d :: Print (ostream & ost) const
{
ost << "SplineCurve: " << points.Size() << " points." << endl;
for (int i = 1; i <= points.Size(); i++)
ost << "P" << i << " = " << points.Get(i) << endl;
}
}
#ifndef FILE_CSG
#define FILE_CSG
/* *************************************************************************/
/* File: geoml.hpp */
/* Author: Joachim Schoeberl */
/* Date: 21. Jun. 98 */
/* *************************************************************************/
#include <myadt.hpp>
#include <gprim.hpp>
#include <meshing.hpp>
namespace netgen
{
#include "surface.hpp"
#include "solid.hpp"
#include "identify.hpp"
#include "singularref.hpp"
#include "csgeom.hpp"
#ifndef SMALLLIB
#define _INCLUDE_MORE
#endif
#ifdef LINUX
#define _INCLUDE_MORE
#endif
#ifdef _INCLUDE_MORE
#include "triapprox.hpp"
#include "algprim.hpp"
#include "brick.hpp"
#include "spline3d.hpp"
#include "manifold.hpp"
#include "curve2d.hpp"
#include "explicitcurve2d.hpp"
#include "gencyl.hpp"
#include "polyhedra.hpp"
#include "extrusion.hpp"
#include "revolution.hpp"
#include "specpoin.hpp"
#include "edgeflw.hpp"
#include "meshsurf.hpp"
#endif
}
#endif
This diff is collapsed.
#ifndef FILE_CSGEOM
#define FILE_CSGEOM
/**************************************************************************/
/* File: csgeom.hh */
/* Author: Joachim Schoeberl */
/* Date: 27. Nov. 97 */
/**************************************************************************/
/**
Constructive Solid Geometry
*/
class TriangleApproximation;
class TATriangle;
/**
A top level object is an entity to be meshed.
I can be either a solid, or one surface patch of a solid.
*/
class TopLevelObject
{
Solid * solid;
Surface * surface;
double red, blue, green;
bool visible, transp;
double maxh;
string material;
int layer;
int bc; // for surface patches, only
public:
TopLevelObject (Solid * asolid,
Surface * asurface = NULL);
const Solid * GetSolid() const { return solid; }
Solid * GetSolid() { return solid; }
const Surface * GetSurface () const { return surface; }
Surface * GetSurface () { return surface; }
void GetData (ostream & ost);
void SetData (istream & ist);
void SetMaxH (double amaxh) { maxh = amaxh; }
double GetMaxH () const { return maxh; }
void SetRGB (double ared, double agreen, double ablue)
{
red = ared;
green = agreen;
blue = ablue;
}
double GetRed () const { return red; }
double GetGreen () const { return green; }
double GetBlue () const { return blue; }
void SetTransparent (bool atransp)
{ transp = atransp; }
bool GetTransparent () const { return transp; }
void SetVisible (bool avisible)
{ visible = avisible; }
bool GetVisible () const { return visible; }
const string GetMaterial () const { return material; }
void SetMaterial (const string & mat) { material = mat; }
int GetLayer () const { return layer; }
void SetLayer (int alayer) { layer = alayer; }
void SetBCProp (int abc) { bc = abc; }
int GetBCProp () const { return bc; }
};
/**
CSGeometry has the whole geometric information
*/
class CSGeometry
{
private:
/// all surfaces
SYMBOLTABLE<Surface*> surfaces;
/// all named solids
SYMBOLTABLE<Solid*> solids;
/// all top level objects: solids and surfaces
ARRAY<TopLevelObject*> toplevelobjects;
/// additional points specified by user
ARRAY<Point<3> > userpoints;
/// triangular approximation of top level objects
ARRAY<TriangleApproximation*> triapprox;
/// increment, if geometry is changed
static int changeval;
/// bounding box of geometry
Box<3> boundingbox;
/// bounding box, if not set by input file
static Box<3> default_boundingbox;
/// identic surfaces are stored by pair of indizes, val = inverse
INDEX_2_HASHTABLE<int> identicsurfaces;
ARRAY<int> isidenticto;
/// identification of boundaries (periodic, thin domains, ...)
/// filename of inputfile
string filename;
public:
CSGeometry ();
CSGeometry (const string & afilename);
~CSGeometry ();
void Clean ();
void Save (ostream & ost);
void Load (istream & ist);
int GetChangeVal() { return changeval; }
void Change() { changeval++; }
void AddSurface (Surface * surf);
void AddSurface (char * name, Surface * surf);
void AddSurfaces (Primitive * prim);
int GetNSurf () const { return surfaces.Size(); }
const Surface * GetSurface (const char * name) const;
const Surface * GetSurface (int i) const
{ return surfaces[i]; }
void SetSolid (const char * name, Solid * sol);
const Solid * GetSolid (const char * name) const;
const Solid * GetSolid (const string & name) const;
int GetNSolids () const { return solids.Size(); }
const Solid * GetSolid (int i) { return solids[i]; }
const SYMBOLTABLE<Solid*> & GetSolids () const { return solids; }
void SetFlags (const char * solidname, const Flags & flags);
int GetNTopLevelObjects () const
{ return toplevelobjects.Size(); }
int SetTopLevelObject (Solid * sol, Surface * surf = NULL);
void GetTopLevelObject (int nr, Solid *& sol, Surface *& surf)
{
sol = toplevelobjects[nr]->GetSolid();
surf = toplevelobjects[nr]->GetSurface();
}
void GetTopLevelObject (int nr, const Solid *& sol, const Surface *& surf) const
{
sol = toplevelobjects[nr]->GetSolid();
surf = toplevelobjects[nr]->GetSurface();
}
TopLevelObject * GetTopLevelObject (const Solid * sol, const Surface * surf = NULL);
TopLevelObject * GetTopLevelObject (int nr)
{ return toplevelobjects[nr]; }
const TopLevelObject * GetTopLevelObject (int nr) const
{ return toplevelobjects[nr]; }
void RemoveTopLevelObject (Solid * sol, Surface * surf = NULL);
void AddUserPoint (const Point<3> & p)
{ userpoints.Append (p); }
int GetNUserPoints () const
{ return userpoints.Size(); }
const Point<3> & GetUserPoint (int nr) const
{ return userpoints[nr]; }
// quick implementations:
ARRAY<SingularFace*> singfaces;
ARRAY<SingularEdge*> singedges;
ARRAY<SingularPoint*> singpoints;
ARRAY<Identification*> identifications;
int GetNIdentifications () { return identifications.Size(); }
void AddIdentification (Identification * ident);
///
void CalcTriangleApproximation(const Box<3> & boundingbox,
double detail, double facets);
///
void FindIdenticSurfaces (double eps);
///
void GetIndependentSurfaceIndices (const Solid * sol,
const BoxSphere<3> & box,
ARRAY<int> & locsurf) const;
///
void GetIndependentSurfaceIndices (const Solid * sol,
const Point<3> & p, Vec<3> & v,
ARRAY<int> & locsurf) const;
///
int GetSurfaceClassRepresentant (int si) const
{ return isidenticto[si]; }
///
const TriangleApproximation * GetTriApprox (int msnr)
{
if (msnr < triapprox.Size())
return triapprox[msnr];
return 0;
}
void IterateAllSolids (SolidIterator & it, int only_once = 0);
void RefineTriangleApprox (Solid * locsol,
int surfind,
const BoxSphere<3> & box,
double detail,
const TATriangle & tria,
TriangleApproximation & tams);
const Box<3> & BoundingBox () const { return boundingbox; }
void SetBoundingBox (const Box<3> & abox)
{
boundingbox = abox;
}
static void SetDefaultBoundingBox (const Box<3> & abox)
{
default_boundingbox = abox;
}
double MaxSize () const;
class BCModification {
public:
int si;
int tlonr;
int bcnr;
};
ARRAY<BCModification> bcmodifications;
};
#endif
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment