mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-27 14:20:34 +05:00
574 lines
11 KiB
C++
574 lines
11 KiB
C++
#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";
|
|
}
|
|
|
|
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 ();
|
|
}
|
|
|
|
|
|
|
|
double Surface :: LocH (const Point<3> & p, double x, double c,
|
|
const MeshingParameters & mparam,
|
|
double hmax) const
|
|
// 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,
|
|
Array<double> & coeffs) const
|
|
{
|
|
classname = "undef";
|
|
coeffs.SetSize (0);
|
|
}
|
|
|
|
void Primitive :: SetPrimitiveData (Array<double> & coeffs)
|
|
{
|
|
;
|
|
}
|
|
|
|
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,
|
|
Array<int> & surfind, double eps) const
|
|
{
|
|
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,
|
|
Array<int> & surfind, double eps) const
|
|
{
|
|
cout << "get tangvecsurfind not implemented" << endl;
|
|
surfind.SetSize (0);
|
|
}
|
|
|
|
void Primitive ::
|
|
GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
|
|
Array<int> & surfind, double eps) const
|
|
{
|
|
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 <= 0)
|
|
return IS_INSIDE;
|
|
else
|
|
return IS_OUTSIDE;
|
|
}
|
|
|
|
|
|
|
|
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;
|
|
}
|
|
}
|
|
|
|
RegisterClassForArchive<Surface> regsurf;
|
|
RegisterClassForArchive<Primitive> regprim;
|
|
RegisterClassForArchive<OneSurfacePrimitive, Surface, Primitive> regosf;
|
|
}
|