ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-04-13 08:37:28 +05:00
Code Issues Packages Projects Releases Wiki Activity

VecInSolid, poly

Browse Source
This commit is contained in:
Joachim Schöberl 2020-10-15 09:29:36 +02:00
parent 54dba89dd8
commit 4f4483794d
2 changed files with 657 additions and 595 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

254
libsrc/csg/polyhedra.cpp
View File

@ -6,10 +6,10 @@
namespace netgen namespace netgen
{ {
Polyhedra::Face::Face (int pi1, int pi2, int pi3, Polyhedra::Face::Face (int pi1, int pi2, int pi3,
const NgArray<Point<3> > & points, const NgArray<Point<3> > & points,
int ainputnr) int ainputnr)
{ {
inputnr = ainputnr; inputnr = ainputnr;
pnums[0] = pi1; pnums[0] = pi1;
@ -43,28 +43,28 @@ Polyhedra::Face::Face (int pi1, int pi2, int pi3,
w1(i) = inv(i,0); w1(i) = inv(i,0);
w2(i) = inv(i,1); w2(i) = inv(i,1);
} }
} }
Polyhedra :: Polyhedra () Polyhedra :: Polyhedra ()
{ {
surfaceactive.SetSize(0); surfaceactive.SetSize(0);
surfaceids.SetSize(0); surfaceids.SetSize(0);
eps_base1 = 1e-8; eps_base1 = 1e-8;
} }
Polyhedra :: ~Polyhedra () Polyhedra :: ~Polyhedra ()
{ {
; ;
} }
Primitive * Polyhedra :: CreateDefault () Primitive * Polyhedra :: CreateDefault ()
{ {
return new Polyhedra(); return new Polyhedra();
} }
INSOLID_TYPE Polyhedra :: BoxInSolid (const BoxSphere<3> & box) const INSOLID_TYPE Polyhedra :: BoxInSolid (const BoxSphere<3> & box) const
{ {
/* /*
for (i = 1; i <= faces.Size(); i++) for (i = 1; i <= faces.Size(); i++)
if (FaceBoxIntersection (i, box)) if (FaceBoxIntersection (i, box))
@ -97,13 +97,13 @@ INSOLID_TYPE Polyhedra :: BoxInSolid (const BoxSphere<3> & box) const
}; };
return PointInSolid (box.Center(), 1e-3 * box.Diam()); return PointInSolid (box.Center(), 1e-3 * box.Diam());
} }
// check how many faces a ray starting in p intersects // check how many faces a ray starting in p intersects
INSOLID_TYPE Polyhedra :: PointInSolid (const Point<3> & p, INSOLID_TYPE Polyhedra :: PointInSolid (const Point<3> & p,
double eps) const double eps) const
{ {
if (!poly_bbox.IsIn (p, eps)) if (!poly_bbox.IsIn (p, eps))
return IS_OUTSIDE; return IS_OUTSIDE;
@ -140,14 +140,14 @@ INSOLID_TYPE Polyhedra :: PointInSolid (const Point<3> & p,
} }
return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE; return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE;
} }
void Polyhedra :: GetTangentialSurfaceIndices (const Point<3> & p, void Polyhedra :: GetTangentialSurfaceIndices (const Point<3> & p,
NgArray<int> & surfind, double eps) const NgArray<int> & surfind, double eps) const
{ {
for (int i = 0; i < faces.Size(); i++) for (int i = 0; i < faces.Size(); i++)
{ {
auto & face = faces[i]; auto & face = faces[i];
@ -166,14 +166,12 @@ void Polyhedra :: GetTangentialSurfaceIndices (const Point<3> & p,
surfind.Append (GetSurfaceId(i)); surfind.Append (GetSurfaceId(i));
} }
} }
#define OLDVECINSOLD INSOLID_TYPE Polyhedra :: VecInSolidOld (const Point<3> & p,
#ifdef OLDVECINSOLD
INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
const Vec<3> & v, const Vec<3> & v,
double eps) const double eps) const
{ {
NgArray<int> point_on_faces; NgArray<int> point_on_faces;
INSOLID_TYPE res(DOES_INTERSECT); INSOLID_TYPE res(DOES_INTERSECT);
@ -238,16 +236,15 @@ INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
// (*testout) << "mindist " << mindist << " res " << res << endl; // (*testout) << "mindist " << mindist << " res " << res << endl;
return res; return res;
} }
#else
// check how many faces a ray starting in p+alpha*v intersects // check how many faces a ray starting in p+alpha*v intersects
INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p, INSOLID_TYPE Polyhedra :: VecInSolidNew (const Point<3> & p,
const Vec<3> & v, const Vec<3> & v,
double eps) const double eps, bool printing) const
{ {
if (!poly_bbox.IsIn (p, eps)) if (!poly_bbox.IsIn (p, eps))
return IS_OUTSIDE; return IS_OUTSIDE;
@ -258,24 +255,35 @@ INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
for (auto & face : faces) for (auto & face : faces)
{ {
Vec<3> v0 = p - points[face.pnums[0]]; Vec<3> v0 = p - points[face.pnums[0]];
if (printing)
{
*testout << "face: ";
for (int j = 0; j < 3; j++)
*testout << points[face.pnums[j]] << " ";
*testout << endl;
}
double lamn = face.nn * v0; double lamn = face.nn * v0;
if (fabs(lamn) < eps) // point is in plance of face if (fabs(lamn) < eps) // point is in plane of face
{ {
double lam1 = face.w1 * v0; double lam1 = face.w1 * v0;
double lam2 = face.w2 * v0; double lam2 = face.w2 * v0;
double lam3 = 1-lam1-lam2; double lam3 = 1-lam1-lam2;
if (printing)
*testout << "lam = " << lam1 << " " << lam2 << " " << lam3 << endl;
if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1) if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1)
{ // point is close to trianlge, perturbe by alpha*v { // point is close to trianlge, perturbe by alpha*v
double dlamn = face.nn*v; double dlamn = face.nn*v;
double dlam1 = face.w1 * v;
double dlam2 = face.w2 * v;
double dlam3 = 1-dlam1-dlam2;
if (fabs(dlamn) < 1e-8) // vec also in plane if (fabs(dlamn) < 1e-8) // vec also in plane
{ {
if (printing)
*testout << "tang in plane" << endl;
double dlam1 = face.w1 * v;
double dlam2 = face.w2 * v;
double dlam3 = 1-dlam1-dlam2;
if (printing)
*testout << "dlam = " << dlam1 << " " << dlam2 << " " << dlam3 << endl;
bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1; bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1;
bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1; bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1;
bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1; bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1;
@ -284,24 +292,41 @@ INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
} }
else // vec out of plane else // vec out of plane
{ {
if (printing)
*testout << "out of plane";
double dlamn = -(face.n * v) / (face.n * n); double dlamn = -(face.n * v) / (face.n * n);
if (printing)
*testout << "dlamn = " << dlamn << endl;
if (dlamn < 0) continue; // ray goes not in direction of face if (dlamn < 0) continue; // ray goes not in direction of face
Vec<3> drs = v0 + lamn * n; Vec<3> drs = v + dlamn * n;
if (printing)
{
*testout << "drs = " << drs << endl;
*testout << "face.w1 = " << face.w1 << endl;
*testout << "face.w2 = " << face.w2 << endl;
}
double dlam1 = dlamn * face.w1 * v; double dlam1 = face.w1 * drs;
double dlam2 = dlamn * face.w2 * v; double dlam2 = face.w2 * drs;
double dlam3 = 1-dlam1-dlam2; double dlam3 = -dlam1-dlam2;
if (printing)
*testout << "dlam = " << dlam1 << " " << dlam2 << " " << dlam3 << endl;
bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1; bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1;
bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1; bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1;
bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1; bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1;
if (in1 && in2 && in3) if (in1 && in2 && in3)
{
if (printing)
*testout << "hit" << endl;
cnt++; cnt++;
} }
} }
} }
}
else else
{ {
double lamn = -(face.n * v0) / (face.n * n); double lamn = -(face.n * v0) / (face.n * n);
@ -314,14 +339,43 @@ INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
double lam2 = face.w2 * rs; double lam2 = face.w2 * rs;
double lam3 = 1-lam1-lam2; double lam3 = 1-lam1-lam2;
if (lam1 >= 0 && lam2 >= 0 && lam3 >= 0) if (lam1 >= 0 && lam2 >= 0 && lam3 >= 0)
{
if (printing)
*testout << "hit" << endl;
cnt++; cnt++;
} }
} }
}
return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE; return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE;
} }
#endif
INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
const Vec<3> & v,
double eps) const
{
return VecInSolidOld (p, v, eps);
/*
auto oldval = VecInSolidOld (p, v, eps);
auto newval = VecInSolidNew (p, v, eps);
if (oldval != newval)
{
*testout << "different decision: oldval = " << oldval
<< " newval = " << newval << endl;
*testout << "p = " << p << ", v = " << v << endl;
VecInSolidNew (p, v, eps, true);
*testout << "Poly:" << endl;
for (auto & face : faces)
{
for (int j = 0; j < 3; j++)
*testout << points[face.pnums[j]] << " ";
*testout << endl;
}
}
return newval;
*/
}
@ -330,12 +384,12 @@ INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
/* /*
INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p, INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
const Vec<3> & v1, const Vec<3> & v1,
const Vec<3> & v2, const Vec<3> & v2,
double eps) const double eps) const
{ {
INSOLID_TYPE res; INSOLID_TYPE res;
res = VecInSolid(p,v1,eps); res = VecInSolid(p,v1,eps);
@ -383,16 +437,16 @@ INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl; cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl;
return Primitive :: VecInSolid2 (p, v1, v2, eps); return Primitive :: VecInSolid2 (p, v1, v2, eps);
} }
*/ */
INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p, INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
const Vec<3> & v1, const Vec<3> & v1,
const Vec<3> & v2, const Vec<3> & v2,
double eps) const double eps) const
{ {
//(*testout) << "VecInSolid2 eps " << eps << endl; //(*testout) << "VecInSolid2 eps " << eps << endl;
INSOLID_TYPE res = VecInSolid(p,v1,eps); INSOLID_TYPE res = VecInSolid(p,v1,eps);
//(*testout) << "VecInSolid = " <<res <<endl; //(*testout) << "VecInSolid = " <<res <<endl;
@ -452,13 +506,13 @@ INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl; cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl;
return Primitive :: VecInSolid2 (p, v1, v2, eps); return Primitive :: VecInSolid2 (p, v1, v2, eps);
} }
void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
NgArray<int> & surfind, double eps) const NgArray<int> & surfind, double eps) const
{ {
Vec<3> v1n = v1; Vec<3> v1n = v1;
v1n.Normalize(); v1n.Normalize();
Vec<3> v2n = v2; // - (v2 * v1n) * v1n; Vec<3> v2n = v2; // - (v2 * v1n) * v1n;
@ -502,7 +556,7 @@ void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec
surfind.Append (GetSurfaceId(faces[i].planenr)); surfind.Append (GetSurfaceId(faces[i].planenr));
} }
} }
} }
@ -515,9 +569,9 @@ void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec
void Polyhedra :: GetPrimitiveData (const char *& classname, void Polyhedra :: GetPrimitiveData (const char *& classname,
NgArray<double> & coeffs) const NgArray<double> & coeffs) const
{ {
classname = "Polyhedra"; classname = "Polyhedra";
coeffs.SetSize(0); coeffs.SetSize(0);
coeffs.Append (points.Size()); coeffs.Append (points.Size());
@ -538,31 +592,31 @@ void Polyhedra :: GetPrimitiveData (const char *& classname,
(*testout) << endl; (*testout) << endl;
} }
*/ */
} }
void Polyhedra :: SetPrimitiveData (NgArray<double> & /* coeffs */) void Polyhedra :: SetPrimitiveData (NgArray<double> & /* coeffs */)
{ {
; ;
} }
void Polyhedra :: Reduce (const BoxSphere<3> & box) void Polyhedra :: Reduce (const BoxSphere<3> & box)
{ {
for (int i = 0; i < planes.Size(); i++) for (int i = 0; i < planes.Size(); i++)
surfaceactive[i] = 0; surfaceactive[i] = 0;
for (int i = 0; i < faces.Size(); i++) for (int i = 0; i < faces.Size(); i++)
if (FaceBoxIntersection (i, box)) if (FaceBoxIntersection (i, box))
surfaceactive[faces[i].planenr] = 1; surfaceactive[faces[i].planenr] = 1;
} }
void Polyhedra :: UnReduce () void Polyhedra :: UnReduce ()
{ {
for (int i = 0; i < planes.Size(); i++) for (int i = 0; i < planes.Size(); i++)
surfaceactive[i] = 1; surfaceactive[i] = 1;
} }
int Polyhedra :: AddPoint (const Point<3> & p) int Polyhedra :: AddPoint (const Point<3> & p)
{ {
if(points.Size() == 0) if(points.Size() == 0)
poly_bbox.Set(p); poly_bbox.Set(p);
else else
@ -570,10 +624,10 @@ int Polyhedra :: AddPoint (const Point<3> & p)
points.Append (p); points.Append (p);
return points.Size(); return points.Size();
} }
int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum) int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum)
{ {
(*testout) << "polyhedra, add face " << pi1 << ", " << pi2 << ", " << pi3 << endl; (*testout) << "polyhedra, add face " << pi1 << ", " << pi2 << ", " << pi3 << endl;
if(pi1 == pi2 || pi2 == pi3 || pi3 == pi1) if(pi1 == pi2 || pi2 == pi3 || pi3 == pi1)
@ -596,20 +650,20 @@ int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum)
n.Normalize(); n.Normalize();
Plane pl (p1, n); Plane pl (p1, n);
// int inverse; // int inverse;
// int identicto = -1; // int identicto = -1;
// for (int i = 0; i < planes.Size(); i++) // for (int i = 0; i < planes.Size(); i++)
// if (pl.IsIdentic (*planes[i], inverse, 1e-9*max3(v1.Length(),v2.Length(),Dist(p2,p3)))) // if (pl.IsIdentic (*planes[i], inverse, 1e-9*max3(v1.Length(),v2.Length(),Dist(p2,p3))))
// { // {
// if (!inverse) // if (!inverse)
// identicto = i; // identicto = i;
// } // }
// // cout << "is identic = " << identicto << endl; // // cout << "is identic = " << identicto << endl;
// identicto = -1; // changed April 10, JS // identicto = -1; // changed April 10, JS
// if (identicto != -1) // if (identicto != -1)
// faces.Last().planenr = identicto; // faces.Last().planenr = identicto;
// else // else
{ {
planes.Append (new Plane (p1, n)); planes.Append (new Plane (p1, n));
surfaceactive.Append (1); surfaceactive.Append (1);
@ -617,15 +671,15 @@ int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum)
faces.Last().planenr = planes.Size()-1; faces.Last().planenr = planes.Size()-1;
} }
// (*testout) << "is plane nr " << faces.Last().planenr << endl; // (*testout) << "is plane nr " << faces.Last().planenr << endl;
return faces.Size(); return faces.Size();
} }
int Polyhedra :: FaceBoxIntersection (int fnr, const BoxSphere<3> & box) const int Polyhedra :: FaceBoxIntersection (int fnr, const BoxSphere<3> & box) const
{ {
/* /*
(*testout) << "check face box intersection, fnr = " << fnr << endl; (*testout) << "check face box intersection, fnr = " << fnr << endl;
(*testout) << "box = " << box << endl; (*testout) << "box = " << box << endl;
@ -656,11 +710,11 @@ int Polyhedra :: FaceBoxIntersection (int fnr, const BoxSphere<3> & box) const
return 1; return 1;
} }
return 0; return 0;
} }
void Polyhedra :: GetPolySurfs(NgArray < NgArray<int> * > & polysurfs) void Polyhedra :: GetPolySurfs(NgArray < NgArray<int> * > & polysurfs)
{ {
int maxnum = -1; int maxnum = -1;
for(int i = 0; i<faces.Size(); i++) for(int i = 0; i<faces.Size(); i++)
@ -675,24 +729,24 @@ void Polyhedra :: GetPolySurfs(NgArray < NgArray<int> * > & polysurfs)
for(int i = 0; i<faces.Size(); i++) for(int i = 0; i<faces.Size(); i++)
polysurfs[faces[i].inputnr]->Append(faces[i].planenr); polysurfs[faces[i].inputnr]->Append(faces[i].planenr);
} }
void Polyhedra::CalcSpecialPoints (NgArray<Point<3> > & pts) const void Polyhedra::CalcSpecialPoints (NgArray<Point<3> > & pts) const
{ {
for (int i = 0; i < points.Size(); i++) for (int i = 0; i < points.Size(); i++)
pts.Append (points[i]); pts.Append (points[i]);
} }
void Polyhedra :: AnalyzeSpecialPoint (const Point<3> & /* pt */, void Polyhedra :: AnalyzeSpecialPoint (const Point<3> & /* pt */,
NgArray<Point<3> > & /* specpts */) const NgArray<Point<3> > & /* specpts */) const
{ {
; ;
} }
Vec<3> Polyhedra :: SpecialPointTangentialVector (const Point<3> & p, int s1, int s2) const Vec<3> Polyhedra :: SpecialPointTangentialVector (const Point<3> & p, int s1, int s2) const
{ {
const double eps = 1e-10*poly_bbox.Diam(); const double eps = 1e-10*poly_bbox.Diam();
for (int fi1 = 0; fi1 < faces.Size(); fi1++) for (int fi1 = 0; fi1 < faces.Size(); fi1++)
@ -800,7 +854,7 @@ Vec<3> Polyhedra :: SpecialPointTangentialVector (const Point<3> & p, int s1, in
} }
return Vec<3> (0,0,0); return Vec<3> (0,0,0);
} }
} }

8
libsrc/csg/polyhedra.hpp
View File

@ -54,10 +54,18 @@ namespace netgen
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const; virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
virtual INSOLID_TYPE PointInSolid (const Point<3> & p, virtual INSOLID_TYPE PointInSolid (const Point<3> & p,
double eps) const; double eps) const;
virtual INSOLID_TYPE VecInSolidNew (const Point<3> & p,
const Vec<3> & v,
double eps, bool printing = false) const;
virtual INSOLID_TYPE VecInSolidOld (const Point<3> & p,
const Vec<3> & v,
double eps) const;
virtual INSOLID_TYPE VecInSolid (const Point<3> & p, virtual INSOLID_TYPE VecInSolid (const Point<3> & p,
const Vec<3> & v, const Vec<3> & v,
double eps) const; double eps) const;
// checks if lim s->0 lim t->0 p + t(v1 + s v2) in solid // checks if lim s->0 lim t->0 p + t(v1 + s v2) in solid
virtual INSOLID_TYPE VecInSolid2 (const Point<3> & p, virtual INSOLID_TYPE VecInSolid2 (const Point<3> & p,
const Vec<3> & v1, const Vec<3> & v1,