diff --git a/libsrc/csg/revolution.cpp b/libsrc/csg/revolution.cpp index 28e96a4f..36f7c614 100644 --- a/libsrc/csg/revolution.cpp +++ b/libsrc/csg/revolution.cpp @@ -948,6 +948,67 @@ namespace netgen return VecInSolid(p,v1+0.01*v2,eps); } + + void Revolution :: + GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, + NgArray & surfind, double eps) const + { + *testout << "tangentialvecsurfind2, p = " << p << endl; + for (int i = 0; i < faces.Size(); i++) + if (faces[i]->PointInFace (p, eps)) + { + *testout << "check face " << i << endl; + Point<2> p2d; + Vec<2> v12d; + faces[i]->CalcProj(p,p2d,v1,v12d); + *testout << "v12d = " << v12d << endl; + auto & spline = faces[i]->GetSpline(); + if (Dist2 (spline.StartPI(), p2d) < sqr(eps)) + { + *testout << "start pi" << endl; + Vec<2> tang = spline.GetTangent(0); + double ip = tang*v12d; + *testout << "ip = " << ip << endl; + if (ip > eps) + surfind.Append(GetSurfaceId(i)); + else if (ip > -eps) + { + Vec<2> v22d; + faces[i]->CalcProj(p,p2d,v2,v22d); + double ip2 = tang*v22d; + *testout << "ip2 = " << ip2 << endl; + if (ip2 > -eps) + surfind.Append(GetSurfaceId(i)); + } + } + else if (Dist2 (faces[i]->GetSpline().EndPI(), p2d) < sqr(eps)) + { + *testout << "end pi" << endl; + + Vec<2> tang = spline.GetTangent(1); + double ip = tang*v12d; + *testout << "ip = " << ip << endl; + if (ip < -eps) + surfind.Append(GetSurfaceId(i)); + else if (ip < eps) + { + Vec<2> v22d; + faces[i]->CalcProj(p,p2d,v2,v22d); + double ip2 = tang*v22d; + *testout << "ip2 = " << ip2 << endl; + if (ip2 < eps) + surfind.Append(GetSurfaceId(i)); + } + } + else + { + *testout << "inner point" << endl; + surfind.Append(GetSurfaceId(i)); + } + } + } + + int Revolution :: GetNSurfaces() const { diff --git a/libsrc/csg/revolution.hpp b/libsrc/csg/revolution.hpp index 8f3d72d5..67f9d9c9 100644 --- a/libsrc/csg/revolution.hpp +++ b/libsrc/csg/revolution.hpp @@ -158,7 +158,10 @@ namespace netgen const Vec<3> & v2, double eps) const; - + virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, + NgArray & surfind, double eps) const; + + virtual int GetNSurfaces() const; virtual Surface & GetSurface (int i = 0); virtual const Surface & GetSurface (int i = 0) const; diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index c3231140..16024e95 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -1908,8 +1908,11 @@ namespace netgen } - if (Abs2 (t) < 1e-8) - continue; + if (Abs2 (t) < 1e-16) + { + cerr << "normal vectors degenerated" << endl; + continue; + } #ifdef DEVELOP *testout << " tangential vector " << t << endl; @@ -1950,10 +1953,15 @@ namespace netgen rhs(0) = -t * (hessej * t); rhs(1) = -t * (hessek * t); - CalcInverse (mat, inv); - t2 = inv * rhs; + CalcInverse (mat, inv); + t2 = inv * rhs; + // t2 = StableSolve (mat, rhs); +#ifdef DEVELOP + *testout << "t = " << t << ", t2 = " << t2 << endl; +#endif - + + /* ageometry.GetIndependentSurfaceIndices (locsol, p, t, surfind2); @@ -2010,7 +2018,9 @@ namespace netgen // locsol2 -> TangentialSolid2 (p, m1, locsol3, surfind3, 1e-9*geomsize); locsol -> TangentialEdgeSolid (p, t, t2, m1, locsol3, surfind3, ideps*geomsize); - +#ifdef DEVELOP + (*testout) << "m1 = " << m1 << ", surfind3 = " << surfind3 << endl; +#endif //ageometry.GetIndependentSurfaceIndices (surfind3); if (surfind3.Contains(surfind2[l])) @@ -2019,6 +2029,9 @@ namespace netgen // locsol2 -> TangentialSolid2 (p, m2, locsol3, surfind3, 1e-9*geomsize); locsol -> TangentialEdgeSolid (p, t, t2, m2, locsol3, surfind3, ideps*geomsize); +#ifdef DEVELOP + (*testout) << "m2 = " << m2 << ", surfind3 = " << surfind3 << endl; +#endif // ageometry.GetIndependentSurfaceIndices (surfind3); diff --git a/libsrc/csg/surface.hpp b/libsrc/csg/surface.hpp index a17fd3da..d2b52cbe 100644 --- a/libsrc/csg/surface.hpp +++ b/libsrc/csg/surface.hpp @@ -293,9 +293,13 @@ namespace netgen const Vec<3> & m, double eps) const; + // for a point p in the surface, into which (closed) surfaces does v point into ? virtual void GetTangentialVecSurfaceIndices (const Point<3> & p, const Vec<3> & v, NgArray & surfind, double eps) const; + // a point p in the surface, and v a tangential vector + // for arbitrary small, but positive t consider q := Project(p+t*v) + // into which (closed) surfaces does v2 point into, when starting from q ? virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, NgArray & surfind, double eps) const; diff --git a/libsrc/gprim/geomfuncs.hpp b/libsrc/gprim/geomfuncs.hpp index ea0070df..c85e0dc6 100644 --- a/libsrc/gprim/geomfuncs.hpp +++ b/libsrc/gprim/geomfuncs.hpp @@ -142,10 +142,24 @@ namespace netgen inline void CalcInverse (const Mat<2,3> & m, Mat<3,2> & inv) { + Vec<3> a0 = m.Row(0); + Vec<3> a1 = m.Row(1); + Vec<3> n = Cross(a0, a1); + Vec<3> d0 = Cross(a1, n); + Vec<3> d1 = Cross(a0, n); + double s0 = 1.0/(a0*d0); + double s1 = 1.0/(a1*d1); + for (int i = 0; i < 3; i++) + { + inv(i,0) = s0*d0(i); + inv(i,1) = s1*d1(i); + } + /* Mat<2,2> a = m * Trans (m); Mat<2,2> ainv; CalcInverse (a, ainv); inv = Trans (m) * ainv; + */ } void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv); @@ -166,6 +180,31 @@ namespace netgen void EigenValues (const Mat<3,3> & m, Vec<3> & ev); void EigenValues (const Mat<2,2> & m, Vec<3> & ev); + + template + Vec<3,T> StableSolve (Mat<2,3,T> mat, Vec<2,T> rhs) + { + Vec<3> a0 = mat.Row(0); + Vec<3> a1 = mat.Row(1); + /* + Vec<3> d = Cross ( Cross (a0, a1), a0); + + double alpha = rhs(0) / a0.Length2(); + double beta = (rhs(1)-alpha* (a0*a1)) / (d*a1); + return alpha * a0 + beta * d; + */ + Vec<3> n = Cross(a0, a1); + Vec<3> d0 = Cross(a1, n); + Vec<3> d1 = Cross(a0, n); + double alpha = rhs(0) / (a0*d0); + double beta = rhs(1) / (a1*d1); + return alpha * d0 + beta * d1; + } + + + + + } #endif