special algorithms in findpoints for degenerated cylinder/plane

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<double> & coeffs) const;
     virtual void SetPrimitiveData (Array<double> & coeffs);
@@ -186,7 +187,9 @@ namespace netgen
   public:
     Cylinder (const Point<3> & aa, const Point<3> & ab, double ar);
     Cylinder (Array<double> & 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<double> & coeffs) const;
     virtual void SetPrimitiveData (Array<double> & coeffs);
     static Primitive * CreateDefault ();

libsrc/csg/python_csg.cpp

@@ -285,6 +285,27 @@ void ExportCSG()
 			self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
           }))
 
+    .def("CloseSurfaces", FunctionPointer
+         ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, int reflevels)
+          {
+            Array<int> 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);
+          }))
+    ;
 }
 
 

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<Plane*>(surf))
+                           {
+                             cnt_plane++;
+                             plane = dynamic_cast<Plane*>(surf);
+                             inv_plane = inv;
+                           }
+                         if (dynamic_cast<Cylinder*>(surf))
+                           {
+                             cnt_cyl++;
+                             cyl = dynamic_cast<Cylinder*>(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<void(Surface*,bool)> & 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<int> & surfind, double eps) const

libsrc/csg/solid.hpp

@@ -73,6 +73,7 @@ namespace netgen
     void GetTangentialSurfaceIndices2 (const Point<3> & p, const Vec<3> & v, Array<int> & surfids, double eps) const;
     void GetTangentialSurfaceIndices3 (const Point<3> & p, const Vec<3> & v, const Vec<3> & v2, Array<int> & surfids, double eps) const;
 
+    void ForEachSurface (const std::function<void(Surface*,bool)> & lambda, bool inv = false) const;
 
     Primitive * GetPrimitive ()
     { return (op == TERM || op == TERM_REF) ? prim : NULL; }

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<int> locsurf;  // attention: array is static
+    // static Array<int> locsurf;  // attention: array is static
+    ArrayMem<int,100> locsurf; 
 
-    static int cntbox = 0;
+    // static int cntbox = 0;
-    cntbox++;
+    // cntbox++;
 
     if (level <= boxesinlevel.Size())
       boxesinlevel.Elem(level)++;
@@ -619,7 +626,9 @@ 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) )