From 01c1411d6546d9daa01acd145ba779be52712c0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Sat, 17 Oct 2020 08:18:06 +0200 Subject: [PATCH] robust implementation of Polyhedra::VecInSolid2 --- libsrc/csg/polyhedra.cpp | 161 ++++++++++++++++++++++++++++++++++++++- libsrc/csg/polyhedra.hpp | 49 +++++++----- 2 files changed, 188 insertions(+), 22 deletions(-) diff --git a/libsrc/csg/polyhedra.cpp b/libsrc/csg/polyhedra.cpp index ae9e53bd..ce35e108 100644 --- a/libsrc/csg/polyhedra.cpp +++ b/libsrc/csg/polyhedra.cpp @@ -355,7 +355,7 @@ namespace netgen const Vec<3> & v, double eps) const { - return VecInSolidOld (p, v, eps); + return VecInSolidNew (p, v, eps); /* auto oldval = VecInSolidOld (p, v, eps); auto newval = VecInSolidNew (p, v, eps); @@ -441,7 +441,8 @@ namespace netgen */ - + // #define OLDVECINSOLID2 +#ifdef OLDVECINSOLID2 INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, @@ -510,6 +511,162 @@ namespace netgen +#else + + + // check how many faces a ray starting in p+alpha*v+alpha^2/2 v2 intersects: + // if p + alpha v is in plane, use v2 + INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p, + const Vec<3> & v, + const Vec<3> & v2, + double eps) const + { + if (!poly_bbox.IsIn (p, eps)) + return IS_OUTSIDE; + + // random (?) direction: + Vec<3> n(-0.424621, 0.1543, 0.89212238); + + int cnt = 0; + for (auto & face : faces) + { + Vec<3> v0 = p - points[face.pnums[0]]; + double lamn = face.nn * v0; + + if (fabs(lamn) < eps) // point is in plane of face + { + double lam1 = face.w1 * v0; + double lam2 = face.w2 * v0; + double lam3 = 1-lam1-lam2; + + if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1) + { // point is close to trianlge, perturbe by alpha*v + double dlamn = face.nn*v; + + if (fabs(dlamn) < 1e-8) // vec also in plane + { + double dlam1 = face.w1 * v; + double dlam2 = face.w2 * v; + double dlam3 = -dlam1-dlam2; + + bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1; + bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1; + bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1; + + // and the same thing for v2 + if (in1 && in2 && in3) + { + double ddlamn = face.nn*v2; + + if (fabs(ddlamn) < 1e-8) // vec2 also in plane + { + double ddlam1 = face.w1 * v2; + double ddlam2 = face.w2 * v2; + double ddlam3 = -ddlam1-ddlam2; + + bool ddin1 = lam1 > eps_base1 || dlam1 > eps_base1 || ddlam1 > -eps_base1; + bool ddin2 = lam2 > eps_base1 || dlam2 > eps_base1 || ddlam2 > -eps_base1; + bool ddin3 = lam3 > eps_base1 || dlam3 > eps_base1 || ddlam3 > -eps_base1; + if (ddin1 && ddin2 && ddin3) + return DOES_INTERSECT; + } + else // vec2 out of plane + { + double ddlamn = -(face.n * v2) / (face.n * n); + if (ddlamn < 0) continue; // ray goes not in direction of face + + Vec<3> drs = v; // + dlamn * n; but dlamn==0 + Vec<3> ddrs = v2 + ddlamn * n; + + double dlam1 = face.w1 * drs; + double dlam2 = face.w2 * drs; + double dlam3 = -dlam1-dlam2; + + double ddlam1 = face.w1 * ddrs; + double ddlam2 = face.w2 * ddrs; + double ddlam3 = -ddlam1-ddlam2; + + bool ddin1 = lam1 > eps_base1 || dlam1 > eps_base1 || ddlam1 > -eps_base1; + bool ddin2 = lam2 > eps_base1 || dlam2 > eps_base1 || ddlam2 > -eps_base1; + bool ddin3 = lam3 > eps_base1 || dlam3 > eps_base1 || ddlam3 > -eps_base1; + + if (ddin1 && ddin2 && ddin3) + cnt++; + } + } + } + else // vec out of plane + { + double dlamn = -(face.n * v) / (face.n * n); + if (dlamn < 0) continue; // ray goes not in direction of face + + Vec<3> drs = v + dlamn * n; + + double dlam1 = face.w1 * drs; + double dlam2 = face.w2 * drs; + double dlam3 = -dlam1-dlam2; + + bool in1 = lam1 > eps_base1 || dlam1 > -eps_base1; + bool in2 = lam2 > eps_base1 || dlam2 > -eps_base1; + bool in3 = lam3 > eps_base1 || dlam3 > -eps_base1; + + if (in1 && in2 && in3) + cnt++; + + } + } + } + else + { + double lamn = -(face.n * v0) / (face.n * n); + + if (lamn < 0) continue; // ray goes not in direction of face + + Vec<3> rs = v0 + lamn * n; + + double lam1 = face.w1 * rs; + double lam2 = face.w2 * rs; + double lam3 = 1-lam1-lam2; + if (lam1 >= 0 && lam2 >= 0 && lam3 >= 0) + cnt++; + } + } + + return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE; + } +#endif + + + + + + + + INSOLID_TYPE Polyhedra :: VecInSolid3 (const Point<3> & p, + const Vec<3> & v1, + const Vec<3> & v2, + double eps) const + { + return VecInSolid2 (p, v1, v2, eps); + } + + INSOLID_TYPE Polyhedra :: VecInSolid4 (const Point<3> & p, + const Vec<3> & v, + const Vec<3> & v2, + const Vec<3> & m, + double eps) const + { + auto res = VecInSolid2 (p, v, v2, eps); + + if (res == DOES_INTERSECT) // following edge second order, let m decide + return VecInSolid2 (p, v, m, eps); + + return res; + } + + + + void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, NgArray & surfind, double eps) const { diff --git a/libsrc/csg/polyhedra.hpp b/libsrc/csg/polyhedra.hpp index 720786e4..4ba183d3 100644 --- a/libsrc/csg/polyhedra.hpp +++ b/libsrc/csg/polyhedra.hpp @@ -48,54 +48,63 @@ namespace netgen public: Polyhedra (); - virtual ~Polyhedra (); + virtual ~Polyhedra () override; static Primitive * CreateDefault (); - virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const; + virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const override; virtual INSOLID_TYPE PointInSolid (const Point<3> & p, - double eps) const; + double eps) const override; 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, const Vec<3> & v, - double eps) const; + double eps) const override; - - - // checks if lim s->0 lim t->0 p + t(v1 + s v2) in solid virtual INSOLID_TYPE VecInSolid2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, - double eps) const; + double eps) const override; + + virtual INSOLID_TYPE VecInSolid3 (const Point<3> & p, + const Vec<3> & v1, + const Vec<3> & v2, + double eps) const override; + virtual INSOLID_TYPE VecInSolid4 (const Point<3> & p, + const Vec<3> & v, + const Vec<3> & v2, + const Vec<3> & m, + double eps) const override; + virtual void GetTangentialSurfaceIndices (const Point<3> & p, - NgArray & surfind, double eps) const; + NgArray & surfind, double eps) const override; virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, - NgArray & surfind, double eps) const; + NgArray & surfind, double eps) const override; - virtual void CalcSpecialPoints (NgArray > & pts) const; + virtual void CalcSpecialPoints (NgArray > & pts) const override; virtual void AnalyzeSpecialPoint (const Point<3> & pt, - NgArray > & specpts) const; - virtual Vec<3> SpecialPointTangentialVector (const Point<3> & p, int s1, int s2) const; + NgArray > & specpts) const override; + virtual Vec<3> SpecialPointTangentialVector (const Point<3> & p, int s1, int s2) const override; - virtual int GetNSurfaces() const + virtual int GetNSurfaces() const override { return planes.Size(); } - virtual Surface & GetSurface (int i) + virtual Surface & GetSurface (int i) override { return *planes[i]; } - virtual const Surface & GetSurface (int i) const + virtual const Surface & GetSurface (int i) const override { return *planes[i]; } - virtual void GetPrimitiveData (const char *& classname, NgArray & coeffs) const; - virtual void SetPrimitiveData (NgArray & coeffs); + virtual void GetPrimitiveData (const char *& classname, NgArray & coeffs) const override; + virtual void SetPrimitiveData (NgArray & coeffs) override; - virtual void Reduce (const BoxSphere<3> & box); - virtual void UnReduce (); + virtual void Reduce (const BoxSphere<3> & box) override; + virtual void UnReduce () override; int AddPoint (const Point<3> & p); int AddFace (int pi1, int pi2, int pi3, int inputnum);