From b81f7f5adaf8ef95189f9f28520ef68839d73b99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Sun, 11 Oct 2020 13:26:34 +0200 Subject: [PATCH] improve robustness for revolution surface meshing (FindSpecialPoints) --- libsrc/csg/revolution.cpp | 22 +++++++- libsrc/csg/revolution.hpp | 5 +- libsrc/csg/specpoin.cpp | 102 +++++++++++++++++++++++++++++++++++--- libsrc/csg/specpoin.hpp | 3 ++ 4 files changed, 123 insertions(+), 9 deletions(-) diff --git a/libsrc/csg/revolution.cpp b/libsrc/csg/revolution.cpp index 57498cb4..28e96a4f 100644 --- a/libsrc/csg/revolution.cpp +++ b/libsrc/csg/revolution.cpp @@ -670,6 +670,23 @@ namespace netgen surfaceactive.Append(1); surfaceids.Append(0); } + + // checking + if (type == 2) + { + auto t0 = spline_in.GetSpline(0).GetTangent(0); + cout << "tstart (must be vertically): " << t0 << endl; + + auto tn = spline_in.GetSpline(nsplines-1).GetTangent(1); + cout << "tend (must be vertically): " << tn << endl; + + for (int i = 0; i < nsplines-1; i++) + { + auto ta = spline_in.GetSpline(i).GetTangent(1); + auto tb = spline_in.GetSpline(i+1).GetTangent(0); + cout << "sin (must not be 0) = " << abs(ta(0)*tb(1)-ta(1)*tb(0)) / (Abs(ta)*Abs(tb)); + } + } } Revolution::~Revolution() @@ -764,8 +781,9 @@ namespace netgen int intersections_before(0), intersections_after(0); double randomx = 7.42357; double randomy = 1.814756; - randomx *= 1./sqrt(randomx*randomx+randomy*randomy); - randomy *= 1./sqrt(randomx*randomx+randomy*randomy); + double randomlen = sqrt(randomx*randomx+randomy*randomy); + randomx *= 1./randomlen; + randomy *= 1./randomlen; const double a = randomy; diff --git a/libsrc/csg/revolution.hpp b/libsrc/csg/revolution.hpp index 0a395e3f..8f3d72d5 100644 --- a/libsrc/csg/revolution.hpp +++ b/libsrc/csg/revolution.hpp @@ -67,7 +67,10 @@ namespace netgen virtual double MaxCurvature () const; //virtual double MaxCurvatureLoc (const Point<3> & /* c */ , // double /* rad */) const; - + + Point<3> P0() const { return p0; } + Vec<3> Axis() const { return v_axis; } + virtual void Project (Point<3> & p) const; virtual Point<3> GetSurfacePoint () const; diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index 1855508c..c3231140 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -463,8 +463,35 @@ namespace netgen return; } - } + if (numprim == 2 && + (typeid(*geometry->GetSurface(locsurf[0])) == typeid(RevolutionFace&)) && + (typeid(*geometry->GetSurface(locsurf[1])) == typeid(RevolutionFace&))) + { + NgArray> pts; + bool check = + ComputeExtremalPoints (static_cast (geometry->GetSurface(locsurf[0])), + static_cast (geometry->GetSurface(locsurf[1])), + pts); + if (check) + { + // int cnt_inbox = 0; + for (auto p : pts) + if (box.IsIn(p)) + { + AddPoint (p, layer); + // cnt_inbox++; + } + return; + /* + if (cnt_inbox == 0) + return; + */ + } + } + + } // end if (numprim <= check_crosspoint) + possiblecrossp = (numprim >= 3) && calccp; @@ -608,9 +635,9 @@ namespace netgen decision = 0; } } - - // (*testout) << "l = " << level << " dec/sureexp = " << decision << sureexp << endl; - +#ifdef DEVELOP + (*testout) << "edgepnt decision = " << decision << " sure = " << sureexp << endl; +#endif if (decision && sureexp) { for (int k1 = 0; k1 < locsurf.Size() - 1; k1++) @@ -890,9 +917,18 @@ namespace netgen f1->CalcGradient (p, g1); f2->CalcGradient (p, g2); - if ( sqr (g1 * g2) > (1 - 1e-10) * Abs2 (g1) * Abs2 (g2)) - return 1; + // if ( sqr (g1 * g2) > (1 - 1e-10) * Abs2 (g1) * Abs2 (g2)) + // return 1; + if ( Abs2 (Cross(g1,g2)) < 1e-10 * Abs2 (g1) * Abs2 (g2)) // same, but stable + { + if (Abs2(vrs) < 1e-12*sqr(size)) // degenerate only if on both surfaces + return 1; + else + return 0; + } + + for (int j = 0; j < 3; j++) { mat(0,j) = g1(j); @@ -1534,7 +1570,61 @@ namespace netgen } + bool SpecialPointCalculation :: + ComputeExtremalPoints (const RevolutionFace * rev1, + const RevolutionFace * rev2, + NgArray > & pts) + { + // if (rev1 -> P0() != rev2 -> P0()) return false; // missing ???? + if (Dist2 (rev1 -> P0(), rev2 -> P0()) > 1e-20*sqr(size)) return false; + // if (rev1 -> Axis() != rev2 -> Axis()) return false; + if ( (rev1 -> Axis()-rev2 -> Axis()).Length2() > 1e-16) return false; + Point<2> p1s = rev1->GetSpline().StartPI(); + Point<2> p1e = rev1->GetSpline().EndPI(); + Point<2> p2s = rev2->GetSpline().StartPI(); + Point<2> p2e = rev2->GetSpline().EndPI(); + + Point<2> p2d; + if (Dist2(p1s,p2e) < 1e-20*sqr(size)) + p2d = p1s; + else if (Dist2(p1e, p2s) < 1e-20*sqr(size)) + p2d = p1e; + else + return false; + *testout << "Norm axis = " << rev1->Axis().Length() << endl; + Point<3> center = rev1->P0() + p2d(0)*rev1->Axis(); + Vec<3> n = rev1->Axis(); + // extremal points of circle, center, normal axis, radius p2d(1) + // Lagrange: + // L(x, lam1, lam2) = x_i + lam1 * (x-c)*v + lam2 * ( |x-c|^2 - r^2 ) + for (double i = 0; i < 3; i++) + { + double lam1 = -n(i) / n.Length2(); + Vec<3> ei(0,0,0); ei(i) = 1; + // double lam2 = 1/(2*p2d(1)) * sqrt(1 - sqr(n(i))/n.Length2()); + double fac = 1-sqr(n(i))/n.Length2(); + // if (fabs(lam2) > 1e-10) + if (fac > 1e-10) + { + double lam2 = 1/(2*p2d(1)) * sqrt(fac); + Point<3> x = center - 1.0/(2*lam2) * (ei + lam1*n); + pts.Append (x); + x = center + 1.0/(2*lam2) * (ei + lam1*n); + pts.Append (x); + + /* + // check: + Point<2> p2d; + rev1 -> CalcProj (x, p2d); + *testout << "special solution, p2d = " << p2d << endl; + rev2 -> CalcProj (x, p2d); + *testout << "special solution, p2d = " << p2d << endl; + */ + } + } + return true; + } diff --git a/libsrc/csg/specpoin.hpp b/libsrc/csg/specpoin.hpp index a4210a15..40e70a82 100644 --- a/libsrc/csg/specpoin.hpp +++ b/libsrc/csg/specpoin.hpp @@ -167,6 +167,9 @@ namespace netgen const Sphere * sphere2, NgArray > & pts); + bool ComputeExtremalPoints (const RevolutionFace * rev1, + const RevolutionFace * rev2, + NgArray > & pts); void ComputeCrossPoints (const Plane * plane1, const Plane * plane2,