netgen/libsrc/csg/surface.cpp

584 lines
11 KiB
C++
Raw Normal View History

2009-01-13 04:40:13 +05:00
#include <mystdlib.h>
#include <myadt.hpp>
#include <csg.hpp>
#include <linalg.hpp>
#include <meshing.hpp>
namespace netgen
{
Surface :: Surface ()
{
maxh = 1e10;
name = new char[7];
strcpy (name, "noname");
bcprop = -1;
bcname = "default";
inverse = false;
2009-01-13 04:40:13 +05:00
}
Surface :: ~Surface()
{
delete [] name;
}
void Surface :: SetName (const char * aname)
{
delete [] name;
name = new char[strlen (aname)+1];
strcpy (name, aname);
}
int Surface :: PointOnSurface (const Point<3> & p,
double eps) const
{
double val = CalcFunctionValue (p);
return fabs (val) < eps;
}
void Surface :: CalcHesse (const Point<3> & point, Mat<3> & hesse) const
{
double dx = 1e-5;
Point<3> hp1, hp2;
Vec<3> g1, g2;
for (int i = 0; i < 3; i++)
{
hp1 = point;
hp2 = point;
hp1(i) += dx;
hp2(i) -= dx;
CalcGradient (hp1, g1);
CalcGradient (hp2, g2);
for (int j = 0; j < 3; j++)
hesse(i, j) = (g1(j) - g2(j)) / (2 * dx);
}
}
/*
void Surface :: GetNormalVector (const Point<3> & p, Vec<3> & n) const
{
CalcGradient (p, n);
n.Normalize();
}
*/
Vec<3> Surface :: GetNormalVector (const Point<3> & p) const
{
Vec<3> n;
CalcGradient (p, n);
n.Normalize();
return n;
}
void Surface :: DefineTangentialPlane (const Point<3> & ap1,
const Point<3> & ap2)
{
p1 = ap1;
p2 = ap2;
ez = GetNormalVector (p1);
ex = p2 - p1;
ex -= (ex * ez) * ez;
ex.Normalize();
ey = Cross (ez, ex);
}
void Surface :: ToPlane (const Point<3> & p3d, Point<2> & pplane,
double h, int & zone) const
{
Vec<3> p1p, n;
n = GetNormalVector (p3d);
if (n * ez < 0)
{
zone = -1;
pplane(0) = 1e8;
pplane(1) = 1e9;
return;
}
p1p = p3d - p1;
pplane(0) = (p1p * ex) / h;
pplane(1) = (p1p * ey) / h;
zone = 0;
}
void Surface :: FromPlane (const Point<2> & pplane,
Point<3> & p3d, double h) const
{
p3d = p1
+ (h * pplane(0)) * ex
+ (h * pplane(1)) * ey;
Project (p3d);
}
void Surface :: Project (Point<3> & p) const
{
Vec<3> n;
double val;
for (int i = 1; i <= 10; i++)
{
val = CalcFunctionValue (p);
if (fabs (val) < 1e-12) return;
CalcGradient (p, n);
p -= (val / Abs2 (n)) * n;
}
}
void Surface :: SkewProject (Point<3> & p, const Vec<3> & direction) const
{
Point<3> startp(p);
double t_old(0),t_new(1);
Vec<3> grad;
for(int i=0; fabs(t_old-t_new) > 1e-20 && i<15; i++)
{
t_old = t_new;
CalcGradient(p,grad);
t_new = t_old - CalcFunctionValue(p)/(grad*direction);
p = startp + t_new*direction;
}
}
double Surface :: MaxCurvature () const
{
return 0.5 * HesseNorm ();
}
double Surface ::
MaxCurvatureLoc (const Point<3> & /* c */ , double /* rad */) const
{
return MaxCurvature ();
}
2014-08-30 06:15:59 +06:00
double Surface :: LocH (const Point<3> & p, double x, double c,
const MeshingParameters & mparam,
double hmax) const
2009-01-13 04:40:13 +05:00
// finds h <= hmax, s.t. h * \kappa_x*h < c
{
/*
double h, hmin, kappa;
hmin = 0;
while (hmin < 0.9 * hmax)
{
h = 0.5 * (hmin + hmax);
kappa = 2 * MaxCurvatureLoc (p, x * h);
if (kappa * h >= c)
hmax = h;
else
hmin = h;
}
return h;
*/
double hret;
double kappa = MaxCurvatureLoc (p, x*hmax);
kappa *= c * mparam.curvaturesafety;
if (hmax * kappa < 1)
hret = hmax;
else
hret = 1 / kappa;
if (maxh < hret)
hret = maxh;
return hret;
}
Primitive :: Primitive ()
{
surfaceids.SetSize (1);
surfaceactive.SetSize (1);
surfaceactive[0] = 1;
}
Primitive :: ~Primitive()
{
;
}
int Primitive :: GetSurfaceId (int i) const
{
return surfaceids[i];
}
void Primitive :: SetSurfaceId (int i, int id)
{
surfaceids[i] = id;
}
void Primitive :: GetPrimitiveData (const char *& classname,
2019-07-09 13:39:16 +05:00
NgArray<double> & coeffs) const
2009-01-13 04:40:13 +05:00
{
classname = "undef";
coeffs.SetSize (0);
}
2019-07-09 13:39:16 +05:00
void Primitive :: SetPrimitiveData (NgArray<double> & coeffs)
2009-01-13 04:40:13 +05:00
{
;
}
Primitive * Primitive :: CreatePrimitive (const char * classname)
{
if (strcmp (classname, "sphere") == 0)
return Sphere::CreateDefault();
if (strcmp (classname, "plane") == 0)
return Plane::CreateDefault();
if (strcmp (classname, "cylinder") == 0)
return Cylinder::CreateDefault();
if (strcmp (classname, "cone") == 0)
return Cone::CreateDefault();
if (strcmp (classname, "brick") == 0)
return Brick::CreateDefault();
stringstream ost;
ost << "Primitve::CreatePrimitive not implemented for " << classname << endl;
throw NgException (ost.str());
}
Primitive * Primitive :: Copy () const
{
stringstream ost;
ost << "Primitve::Copy not implemented for " << typeid(*this).name() << endl;
throw NgException (ost.str());
}
void Primitive :: Transform (Transformation<3> & trans)
{
stringstream ost;
ost << "Primitve::Transform not implemented for " << typeid(*this).name() << endl;
throw NgException (ost.str());
}
void Primitive :: GetTangentialSurfaceIndices (const Point<3> & p,
2019-07-09 13:39:16 +05:00
NgArray<int> & surfind, double eps) const
2009-01-13 04:40:13 +05:00
{
for (int j = 0; j < GetNSurfaces(); j++)
if (fabs (GetSurface(j).CalcFunctionValue (p)) < eps)
if (!surfind.Contains (GetSurfaceId(j)))
surfind.Append (GetSurfaceId(j));
}
void Primitive ::
GetTangentialVecSurfaceIndices (const Point<3> & p, const Vec<3> & v,
2019-07-09 13:39:16 +05:00
NgArray<int> & surfind, double eps) const
2009-01-13 04:40:13 +05:00
{
cout << "get tangvecsurfind not implemented" << endl;
surfind.SetSize (0);
}
void Primitive ::
GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
2019-07-09 13:39:16 +05:00
NgArray<int> & surfind, double eps) const
2009-01-13 04:40:13 +05:00
{
for (int j = 0; j < GetNSurfaces(); j++)
{
if (fabs (GetSurface(j).CalcFunctionValue (p)) < eps)
{
Vec<3> grad;
GetSurface(j).CalcGradient (p, grad);
if (sqr (grad * v1) < 1e-6 * v1.Length2() * grad.Length2() &&
sqr (grad * v2) < 1e-6 * v2.Length2() * grad.Length2() ) // new, 18032006 JS
{
if (!surfind.Contains (GetSurfaceId(j)))
surfind.Append (GetSurfaceId(j));
}
}
}
}
INSOLID_TYPE Primitive ::
VecInSolid2 (const Point<3> & p,
const Vec<3> & v1,
const Vec<3> & v2,
double eps) const
{
//(*testout) << "Primitive::VecInSolid2" << endl;
Point<3> hp = p + 1e-3 * v1 + 1e-5 * v2;
INSOLID_TYPE res = PointInSolid (hp, eps);
// (*testout) << "vectorin2, type = " << typeid(*this).name() << ", res = " << res << endl;
return res;
}
INSOLID_TYPE Primitive ::
VecInSolid3 (const Point<3> & p,
const Vec<3> & v1,
const Vec<3> & v2,
double eps) const
{
//(*testout) << "Primitive::VecInSolid3" << endl;
return VecInSolid (p, v1, eps);
}
INSOLID_TYPE Primitive ::
VecInSolid4 (const Point<3> & p,
const Vec<3> & v,
const Vec<3> & v2,
const Vec<3> & m,
double eps) const
{
return VecInSolid2 (p, v, m, eps);
}
OneSurfacePrimitive :: OneSurfacePrimitive()
{
;
}
OneSurfacePrimitive :: ~OneSurfacePrimitive()
{
;
}
INSOLID_TYPE OneSurfacePrimitive ::
PointInSolid (const Point<3> & p,
double eps) const
{
double hv1 = (GetSurface(0).CalcFunctionValue(p));
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
return DOES_INTERSECT;
}
INSOLID_TYPE OneSurfacePrimitive ::
VecInSolid (const Point<3> & p, const Vec<3> & v,
double eps) const
{
double hv1 = (GetSurface(0).CalcFunctionValue(p));
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
Vec<3> hv;
GetSurface(0).CalcGradient (p, hv);
hv1 = v * hv;
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
return DOES_INTERSECT;
}
INSOLID_TYPE OneSurfacePrimitive ::
VecInSolid2 (const Point<3> & p,
const Vec<3> & v1,
const Vec<3> & v2,
double eps) const
{
double hv1 = (GetSurface(0).CalcFunctionValue(p));
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
Vec<3> hv;
GetSurface(0).CalcGradient (p, hv);
hv1 = v1 * hv;
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
double hv2 = v2 * hv;
if (hv2 <= -eps)
return IS_INSIDE;
if (hv2 >= eps)
return IS_OUTSIDE;
return DOES_INTERSECT;
/*
2009-01-13 04:40:13 +05:00
double hv2 = v2 * hv;
if (hv2 <= 0)
return IS_INSIDE;
else
return IS_OUTSIDE;
*/
2009-01-13 04:40:13 +05:00
}
INSOLID_TYPE OneSurfacePrimitive ::
VecInSolid3 (const Point<3> & p, const Vec<3> & v, const Vec<3> & v2,
double eps) const
{
//(*testout) << "OneSurfacePrimitive::VecInSolid3" << endl;
double hv1 = (GetSurface(0).CalcFunctionValue(p));
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
Vec<3> grad;
GetSurface(0).CalcGradient (p, grad);
hv1 = v * grad;
if (hv1 <= -eps) return IS_INSIDE;
if (hv1 >= eps) return IS_OUTSIDE;
Mat<3> hesse;
GetSurface(0).CalcHesse (p, hesse);
double hv2 = v2 * grad + v * (hesse * v);
if (hv2 <= -eps) return IS_INSIDE;
if (hv2 >= eps) return IS_OUTSIDE;
return DOES_INTERSECT;
}
INSOLID_TYPE OneSurfacePrimitive ::
VecInSolid4 (const Point<3> & p, const Vec<3> & v, const Vec<3> & v2,
const Vec<3> & m,
double eps) const
{
double hv1 = (GetSurface(0).CalcFunctionValue(p));
if (hv1 <= -eps)
return IS_INSIDE;
if (hv1 >= eps)
return IS_OUTSIDE;
Vec<3> grad;
GetSurface(0).CalcGradient (p, grad);
hv1 = v * grad;
if (hv1 <= -eps) return IS_INSIDE;
if (hv1 >= eps) return IS_OUTSIDE;
Mat<3> hesse;
GetSurface(0).CalcHesse (p, hesse);
double hv2 = v2 * grad + v * (hesse * v);
if (hv2 <= -eps) return IS_INSIDE;
if (hv2 >= eps) return IS_OUTSIDE;
double hv3 = m * grad;
if (hv3 <= -eps) return IS_INSIDE;
if (hv3 >= eps) return IS_OUTSIDE;
return DOES_INTERSECT;
}
int OneSurfacePrimitive :: GetNSurfaces() const
{
return 1;
}
Surface & OneSurfacePrimitive :: GetSurface (int i)
{
return *this;
}
const Surface & OneSurfacePrimitive :: GetSurface (int i) const
{
return *this;
}
void ProjectToEdge (const Surface * f1, const Surface * f2, Point<3> & hp)
{
Vec<2> rs, lam;
Vec<3> a1, a2;
Mat<2> a;
int i = 10;
while (i > 0)
{
i--;
rs(0) = f1 -> CalcFunctionValue (hp);
rs(1) = f2 -> CalcFunctionValue (hp);
f1->CalcGradient (hp, a1);
f2->CalcGradient (hp, a2);
double alpha = fabs(a1*a2)/sqrt(a1.Length2()*a2.Length2());
if(fabs(1.-alpha) < 1e-6)
{
if(fabs(rs(0)) >= fabs(rs(1)))
f1 -> Project(hp);
else
f2 -> Project(hp);
}
else
{
a(0,0) = a1 * a1;
a(0,1) = a(1,0) = a1 * a2;
a(1,1) = a2 * a2;
a.Solve (rs, lam);
hp -= lam(0) * a1 + lam(1) * a2;
}
if (Abs2 (rs) < 1e-24 && i > 1) i = 1;
}
}
2018-12-06 21:53:44 +05:00
RegisterClassForArchive<Surface> regsurf;
RegisterClassForArchive<Primitive> regprim;
RegisterClassForArchive<OneSurfacePrimitive, Surface, Primitive> regosf;
2009-01-13 04:40:13 +05:00
}