From 2350e8b9d55b2434a878b4e2bf898641c924e4d5 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 7 Dec 2014 09:28:39 +0000 Subject: [PATCH] special algorithms in findpoints for degenerated cylinder/plane --- libsrc/csg/algprim.hpp | 7 ++- libsrc/csg/python_csg.cpp | 36 ++++++++++-- libsrc/csg/solid.cpp | 112 ++++++++++++++++++++++++++++++++++++++ libsrc/csg/solid.hpp | 1 + libsrc/csg/specpoin.cpp | 19 +++++-- 5 files changed, 164 insertions(+), 11 deletions(-) diff --git a/libsrc/csg/algprim.hpp b/libsrc/csg/algprim.hpp index c71dc8c2..0459108a 100644 --- a/libsrc/csg/algprim.hpp +++ b/libsrc/csg/algprim.hpp @@ -63,7 +63,8 @@ namespace netgen public: /// Plane (const Point<3> & ap, Vec<3> an); - + Point<3> P() const { return p; } + Vec<3> N() const { return n; } virtual void GetPrimitiveData (const char *& classname, Array & coeffs) const; virtual void SetPrimitiveData (Array & coeffs); @@ -186,7 +187,9 @@ namespace netgen public: Cylinder (const Point<3> & aa, const Point<3> & ab, double ar); Cylinder (Array & coeffs); - + Point<3> A() const { return a; } + Point<3> B() const { return b; } + double R() const { return r; } virtual void GetPrimitiveData (const char *& classname, Array & coeffs) const; virtual void SetPrimitiveData (Array & coeffs); static Primitive * CreateDefault (); diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 86658d2c..145c65fa 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -111,9 +111,9 @@ public: void SetColor(double ared, double agreen, double ablue) { - red = ared; - green = agreen; - blue = ablue; + red = ared; + green = agreen; + blue = ablue; } double GetRed() const { return red; } @@ -285,6 +285,27 @@ void ExportCSG() self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent()); })) + .def("CloseSurfaces", FunctionPointer + ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, int reflevels) + { + Array si1, si2; + s1->GetSolid()->GetSurfaceIndices (si1); + s2->GetSolid()->GetSurfaceIndices (si2); + cout << "surface ids1 = " << si1 << endl; + cout << "surface ids2 = " << si2 << endl; + + Flags flags; + const TopLevelObject * domain = nullptr; + self.AddIdentification + (new CloseSurfaceIdentification + (self.GetNIdentifications()+1, self, + self.GetSurface (si1[0]), self.GetSurface (si2[0]), + domain, + flags)); + }), + (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"), bp::arg("reflevels")=2) + ) + .add_property ("ntlo", &CSGeometry::GetNTopLevelObjects) ; @@ -300,7 +321,14 @@ void ExportCSG() return dummy; })) ; - + bp::def("ZRefinement", FunctionPointer + ([](Mesh & mesh, CSGeometry & geom) + { + ZRefinementOptions opt; + opt.minref = 5; + ZRefinement (mesh, &geom, opt); + })) + ; } diff --git a/libsrc/csg/solid.cpp b/libsrc/csg/solid.cpp index aba39d8a..4193db88 100644 --- a/libsrc/csg/solid.cpp +++ b/libsrc/csg/solid.cpp @@ -1231,6 +1231,88 @@ namespace netgen Solid * Solid :: RecGetReducedSolid (const BoxSphere<3> & box, INSOLID_TYPE & in) const { + + + { + // checking special case for degenerated plane - cylinder, Dec 2014 + int cnt_plane = 0, cnt_cyl = 0; + bool inv_plane, inv_cyl; + Plane * plane; + Cylinder * cyl; + + ForEachSurface ( [&] (Surface * surf, bool inv) + { + if (dynamic_cast(surf)) + { + cnt_plane++; + plane = dynamic_cast(surf); + inv_plane = inv; + } + if (dynamic_cast(surf)) + { + cnt_cyl++; + cyl = dynamic_cast(surf); + inv_cyl = inv; + } + }); + + if (cnt_plane == 1 && cnt_cyl == 1) + { + double scala = (cyl->A()-plane->P()) * plane->N(); + double scalb = (cyl->B()-plane->P()) * plane->N(); + double scal = plane->N() * plane->N(); + if ( ( fabs (scala*scala - cyl->R()*cyl->R()*scal) < 1e-10*cyl->R()*cyl->R() ) && + ( fabs (scalb*scalb - cyl->R()*cyl->R()*scal) < 1e-10*cyl->R()*cyl->R() ) ) + { + // intersection edge in box ? + Point<3> p0 = cyl->A() - (scala/scal) * plane->N(); + Vec<3> vedge = cyl->B() - cyl->A(); + Vec<3> ve_center = box.Center()-p0; + + // dist(lam) = \| ve_center \|^2 - 2 lam (vedge, ve_center) + lam^2 \| vedge \|^2 + + double num = vedge*ve_center; + double den = vedge*vedge; + + double dist_edge_center2 = ve_center*ve_center - num * num /den; + + + bool edge_in_box = dist_edge_center2 < sqr (box.Diam()); + + + if (!edge_in_box) + { + if (op == SECTION) + { + // cout << "solid = " << *this << endl; + if (!inv_cyl && !inv_plane && scala < 0) + { + cout << "fix for degenerated cyl-plane edge: just the cylinder" << endl; + Solid * sol = new Solid (cyl); + sol -> op = TERM_REF; + return sol; + } + } + + if (op == UNION) + { + // cout << "solid = " << *this << ", inv_plane = " << inv_plane << " inv_cyl = " << inv_cyl << " scalb " << scalb << endl; + if (!inv_plane && !inv_cyl && (scala < 0)) + { + cout << "fix for degenerated cyl-plane edge: just the plane" << endl; + // return new Solid (plane); + Solid * sol = new Solid (plane); + sol -> op = TERM_REF; + return sol; + } + } + ; // *testout << "have 1 plane and 1 cyl, degenerated" << endl; + } + } + } + } + + Solid * redsol = NULL; switch (op) @@ -1417,6 +1499,36 @@ namespace netgen } } + void Solid :: ForEachSurface (const std::function & lambda, bool inv) const + { + switch (op) + { + case TERM: case TERM_REF: + { + for (int j = 0; j < prim->GetNSurfaces(); j++) + if (prim->SurfaceActive (j)) + lambda (&prim->GetSurface(j), inv); + break; + } + case UNION: + case SECTION: + { + s1 -> ForEachSurface (lambda, inv); + s2 -> ForEachSurface (lambda, inv); + break; + } + case SUB: + { + s1 -> ForEachSurface (lambda, !inv); + break; + } + case ROOT: + { + s1 -> ForEachSurface (lambda, inv); + break; + } + } + } void Solid :: GetTangentialSurfaceIndices (const Point<3> & p, Array & surfind, double eps) const diff --git a/libsrc/csg/solid.hpp b/libsrc/csg/solid.hpp index ebd7c88c..17e03512 100644 --- a/libsrc/csg/solid.hpp +++ b/libsrc/csg/solid.hpp @@ -73,6 +73,7 @@ namespace netgen void GetTangentialSurfaceIndices2 (const Point<3> & p, const Vec<3> & v, Array & surfids, double eps) const; void GetTangentialSurfaceIndices3 (const Point<3> & p, const Vec<3> & v, const Vec<3> & v2, Array & surfids, double eps) const; + void ForEachSurface (const std::function & lambda, bool inv = false) const; Primitive * GetPrimitive () { return (op == TERM || op == TERM_REF) ? prim : NULL; } diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index f4c5f946..10da4bb9 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -199,15 +199,22 @@ namespace netgen throw NgException (err.c_str()); } + if (level == 40 || level == 41 || level == 45) + { + *testout << "level = " << level << " cp = " << calccp << " ep = " << calcep << ", box = " << box << ", solid = " << *sol << endl; + } + + bool decision; bool possiblecrossp, possibleexp; // possible cross or extremalpoint bool surecrossp = 0, sureexp = 0; // sure ... - static Array locsurf; // attention: array is static + // static Array locsurf; // attention: array is static + ArrayMem locsurf; - static int cntbox = 0; - cntbox++; + // static int cntbox = 0; + // cntbox++; if (level <= boxesinlevel.Size()) boxesinlevel.Elem(level)++; @@ -619,8 +626,10 @@ namespace netgen Point<3> pp; if (IsEdgeExtremalPoint (surf1, surf2, p, pp, box.Diam()/2)) { - (*testout) << "extremalpoint (nearly) found:" << pp << endl; - + (*testout) << "extremalpoint (nearly) found:" << pp + << "box.diam = " << box.Diam() << ", dist = " << Dist(pp,box.Center()) + << endl; + if (Dist (pp, box.Center()) < box.Diam()/2 && sol -> IsIn (pp, 1e-6*size) && !sol->IsStrictIn (pp, 1e-6*size) ) {