further fixes

libsrc/csg/edgeflw.cpp

@@ -44,6 +44,7 @@ namespace netgen
     for (PointIndex pi : mesh.Points().Range())    
       meshpoint_tree->Insert (mesh[pi], pi);
 
+
     // add all special points before edge points (important for periodic identification)
     // JS, Jan 2007
     const double di=1e-7*geometry.MaxSize();
@@ -454,7 +455,6 @@ namespace netgen
 		refedges[i].surfnr2 = geometry.GetSurfaceClassRepresentant(refedges[i].surfnr2);
 	      }
 
-
 	    /*
 	      for (int i = oldnseg+1; i <= mesh.GetNSeg(); i++)
 	      for (int j = 1; j <= oldnseg; j++)
@@ -485,7 +485,6 @@ namespace netgen
 		      layer,
 		      mesh);
 	  }
-
 	for(int i=0; i<refedges.Size(); i++)
 	  {
 	    auto splinesurface = dynamic_cast<const SplineSurface*>(geometry.GetSurface(refedges[i].surfnr1));
@@ -494,6 +493,17 @@ namespace netgen
 		auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
 		mesh.SetCD2Name(refedges[i].edgenr,*name);
 	      }
+	    else
+	      {
+		auto splinesurface2 = dynamic_cast<const SplineSurface*>(geometry.GetSurface(refedges[i].surfnr2));
+	    if(splinesurface2)
+	      {
+		auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
+		mesh.SetCD2Name(refedges[i].edgenr,*name);
+	      }
+		
+	      }
+	    
 	  }
 
 	/*

libsrc/csg/python_csg.cpp

@@ -252,7 +252,28 @@ DLL_HEADER void ExportCSG()
            }))
     ;
 
-  bp::class_<SplineSurface, shared_ptr<SplineSurface>,boost::noncopyable> ("SplineSurface")
+  bp::class_<SplineSurface, shared_ptr<SplineSurface>,boost::noncopyable> ("SplineSurface",
+			 "A surface for co dim 2 integrals on the splines", bp::no_init)
+    .def("__init__", bp::make_constructor (FunctionPointer
+					   ([](const shared_ptr<SPSolid> base, bp::list cuts)
+	     {
+	       auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive());
+	       auto acuts = new Array<OneSurfacePrimitive*>();
+	       for(int i = 0; i<bp::len(cuts);i++)
+		 {
+		   bp::extract<shared_ptr<SPSolid>> sps(cuts[i]);
+		   if(!sps.check())
+		     throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
+		   auto sp = dynamic_cast<OneSurfacePrimitive*>(sps()->GetSolid()->GetPrimitive());
+		   if(sp)
+		     acuts->Append(sp);
+		   else
+		     throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
+		 }
+	       if(!primitive)
+		 throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!");
+	       return make_shared<SplineSurface>(primitive,acuts);
+	     }),bp::default_call_policies(),(bp::arg("base"), bp::arg("cuts")=bp::list())))
     .def("AddPoint", FunctionPointer
 	 ([] (SplineSurface & self, double x, double y, double z, bool hpref)
 	  {
@@ -261,12 +282,12 @@ DLL_HEADER void ExportCSG()
 	  }),
 	 (bp::arg("self"),bp::arg("x"),bp::arg("y"),bp::arg("z"),bp::arg("hpref")=false))
     .def("AddSegment", FunctionPointer
-	 ([] (SplineSurface & self, int i1, int i2, string bcname)
+	 ([] (SplineSurface & self, int i1, int i2, string bcname, double maxh)
 	  {
 	    auto str = new string(bcname);
-	    self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str);
+	    self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str,maxh);
 	  }),
-	 (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default"))
+	 (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default", bp::arg("maxh")=-1.))
     ;
   
 #if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100)
@@ -450,13 +471,15 @@ DLL_HEADER void ExportCSG()
     .def("AddSplineSurface", FunctionPointer
 	 ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
 	  {
-	    auto planes = surf->CreatePlanes();
+	    auto cuttings = surf->CreateCuttingSurfaces();
 	    auto spsol = make_shared<SPSolid>(new Solid(&*surf));
-	    for(auto plane : (*planes)){
+	    for(auto cut : (*cuttings)){
-	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(plane)));
+	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut)));
 	    }
 	    spsol->AddSurfaces(self);
 	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf);
+	    for(auto p : surf->GetPoints())
+		self.AddUserPoint(p);
 	  }),
 	 (bp::arg("self"), bp::arg("SplineSurface")))
 

libsrc/csg/splinesurface.cpp

@@ -3,47 +3,60 @@
 
 namespace netgen
 {
-  SplineSurface :: SplineSurface() : OneSurfacePrimitive()
-{ ; }
-
-SplineSurface :: ~SplineSurface() { ; }
-
 void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref)
 {
-  geompoints.Append(GeomPoint<3>(p,reffac));
+  auto pp = Point<3>(p);
+  geompoints.Append(GeomPoint<3>(pp,reffac));
   geompoints.Last().hpref = hpref;
 }
   
-  void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname)
+  void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh)
   {
-      splines.Append(spline);
+    splines.Append(spline);
-      bcnames.Append(bcname);
+    bcnames.Append(bcname);
+    maxh.Append(amaxh);
   }
 
   string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const
   {
-    (*testout) << "segment: " << p1 << ", " << p2 << endl;
+    
+    double eps = 1e-5;
     for(int i=0; i<splines.Size(); i++)
       {
-	(*testout) << "spline: " << splines[i]->GetPoint(0) << ", " << splines[i]->GetPoint(1) << endl;
+	auto pp1 = Point<3>(splines[i]->GetPoint(0));
-	if (((splines[i]->GetPoint(0)-p1).Length()<1e-12 && (splines[i]->GetPoint(1)-p2).Length() < 1e-12) || ((splines[i]->GetPoint(0)-p2).Length() < 1e-12 && (splines[i]->GetPoint(1)-p1).Length() < 1e-12))
+	Project(pp1);
+	auto pp2 = Point<3>(splines[i]->GetPoint(1));
+	Project(pp2);
+	if (((pp1-p1).Length()<eps && (pp2-p2).Length() < eps) || ((pp1-p2).Length() < eps && (pp2-p1).Length() < eps))
 	  {
-	    (*testout) << "return bcname: " << *bcnames[i] << endl;
 	    return bcnames[i];
 	  }
       }
     return new string("default");
   }
 
-  Array<Plane*>* SplineSurface :: CreatePlanes() const
+  Array<OneSurfacePrimitive*>* SplineSurface :: CreateCuttingSurfaces() const
   {
-    auto planes = new Array<Plane*>();
+    auto cuttings = new Array<OneSurfacePrimitive*>();
-    auto sol = new Solid(new Plane(splines[0]->GetPoint(0),Cross(splines[0]->GetTangent(0),GetNormalVector(splines[0]->GetPoint(0)))));
+    for (auto cut : *cuts)
-    for(auto spline : splines)
+      cuttings->Append(cut);
+    for(int i = 0; i<splines.Size(); i++)
       {
-	planes->Append(new Plane(spline->GetPoint(0),-spline->GetTangent(0)));
+	auto spline = splines[i];
+	auto lineseg = dynamic_cast<LineSeg<3>*>(spline);
+	auto p1 = Point<3>(spline->GetPoint(0));
+	Project(p1);
+	auto p2 = Point<3>(spline->GetPoint(1));
+	Project(p2);
+	auto vec = Vec<3>(p2)-Vec<3>(p1);
+	auto plane = new Plane(p1,-Cross(vec,baseprimitive->GetNormalVector(p1)));
+	if(maxh[i]>0)
+	  {
+	  plane->SetMaxH(maxh[i]);
+	  }
+	cuttings->Append(plane);
       }
-    return planes;
+    return cuttings;
   }
   
   void SplineSurface :: Print(ostream & str) const
@@ -51,89 +64,4 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const
   str << "SplineSurface " << endl;
 }
 
-  void SplineSurface :: Project(Point<3> & p3d) const
-  {
-    double val = CalcFunctionValue(p3d);
-    p3d -= val * GetNormalVector(p3d);
-  }
-
-  
-double SplineSurface :: CalcFunctionValue (const Point<3> & point) const
-{
-  auto v1 = splines[0]->GetTangent(0);
-  auto v2 = splines.Last()->GetTangent(1);
-  auto p1 = splines[0]->GetPoint(0);
-  auto n = Cross(v1,v2);
-  n.Normalize();
-  return n * (point-p1);
-}
-
-void SplineSurface :: CalcGradient (const Point<3> & point, Vec<3> & grad) const
-{
-  auto v1 = splines[0]->GetTangent(0);
-  auto v2 = splines.Last()->GetTangent(1);
-  auto p1 = splines[0]->GetPoint(0);
-  grad = Cross(v1,v2);
-  grad.Normalize();
-}
-
-  double SplineSurface :: HesseNorm() const
-{
-  return 0;
-}
-
-  Point<3> SplineSurface :: GetSurfacePoint() const
-{
-  return splines[0]->GetPoint(0);
-}
-
-  
-  void SplineSurface :: CalcSpecialPoints(Array<Point<3>> & pts) const
-  {
-    for(auto pt : geompoints)
-      {
-	pts.Append(Point<3>(pt));
-      }
-  }
-
-
-
-INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const
-{
-  
-  int i;
-  double val;
-   
-  val = CalcFunctionValue (box.Center());
-  if (val > box.Diam() / 2) return IS_OUTSIDE;
-  if (val < -box.Diam() / 2) return IS_INSIDE;
-  
-  Vec<3> n = GetNormalVector(box.Center());
-  double cx = n[0];
-  double cy = n[1];
-  double cz = n[2];
-  
-  if (val > 0)
-    {
-      Vec<3> vdiag = box.PMax() - box.PMin();
-      double modify = (vdiag(0) * fabs (cx) + 
-		       vdiag(1) * fabs (cy) + 
-		       vdiag(2) * fabs (cz) ) / 2;
-      
-      if (val - modify < 0)
-	return DOES_INTERSECT;
-      return IS_OUTSIDE;
-    }
-  else
-    {
-      Vec<3> vdiag = box.PMax() - box.PMin();
-      double modify =  (vdiag(0) * fabs (cx) + 
-			vdiag(1) * fabs (cy) + 
-			vdiag(2) * fabs (cz) ) / 2;
-      if (val + modify > 0)
-	return DOES_INTERSECT;
-      return IS_INSIDE;
-    }
-}
-  
 }

libsrc/csg/splinesurface.hpp

@@ -4,19 +4,25 @@
 
 namespace netgen
 {
-  class SplineSurface : public  OneSurfacePrimitive
+  class SplineSurface  : public  OneSurfacePrimitive
   {
   protected:
     Array<GeomPoint<3>> geompoints;
     Array<SplineSeg<3>*> splines;
     Array<string*> bcnames;
+    Array<double> maxh;
+    OneSurfacePrimitive* baseprimitive;
+    Array<OneSurfacePrimitive*>* cuts;
     
   public:
-    SplineSurface();
+    SplineSurface(OneSurfacePrimitive* abaseprimitive, Array<OneSurfacePrimitive*>* acuts) :
-    virtual ~SplineSurface();
+      OneSurfacePrimitive(), baseprimitive(abaseprimitive), cuts(acuts)
+    { ; }
+    virtual ~SplineSurface() { ; }
     
     const Array<SplineSeg<3>*> & GetSplines() const { return splines; }
     int GetNSplines() const { return splines.Size(); }
+    const Array<GeomPoint<3>>& GetPoints() const { return geompoints; }
     string GetSplineType(const int i) const { return splines[i]->GetType(); }
     SplineSeg<3> & GetSpline(const int i) { return *splines[i]; }
     const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; }
@@ -26,9 +32,30 @@ namespace netgen
     string* GetBCNameOf(Point<3> p1, Point<3> p2) const;
     
     DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false);
-    void AppendSegment(SplineSeg<3>* spline, string* bcname);
+    void AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh = -1);
-    Array<Plane*>* CreatePlanes() const;
 
+    OneSurfacePrimitive* GetBase() const { return baseprimitive; } 
+
+    Array<OneSurfacePrimitive*>* CreateCuttingSurfaces() const;
+
+    
+    virtual void Project (Point<3> & p3d) const { baseprimitive->Project(p3d); }
+    virtual double CalcFunctionValue (const Point<3> & point) const
+    { return baseprimitive->CalcFunctionValue (point); }
+    virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const
+    { baseprimitive->CalcGradient (point,grad); }
+    virtual double HesseNorm () const
+    { return baseprimitive->HesseNorm(); }
+    virtual Point<3> GetSurfacePoint () const
+    { return baseprimitive->GetSurfacePoint(); }
+    virtual void CalcSpecialPoints(Array<Point<3>> & pts) const
+    { baseprimitive->CalcSpecialPoints(pts); }
+
+    virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const
+    { return baseprimitive->BoxInSolid(box); }
+    
+    
+    /*
     virtual void Project (Point<3> & p3d) const;
     virtual double CalcFunctionValue (const Point<3> & point) const;
     virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const;
@@ -38,13 +65,15 @@ namespace netgen
 
     virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const;
 
-    
+    */
     
     virtual void Print (ostream & str) const;
     
     
   };
 
+
+
 }
 
 #endif