first push

libsrc/csg/CMakeLists.txt

@@ -5,7 +5,7 @@ add_library(csg ${NG_LIB_TYPE}
         explicitcurve2d.cpp extrusion.cpp gencyl.cpp genmesh.cpp identify.cpp  
         manifold.cpp meshsurf.cpp polyhedra.cpp revolution.cpp singularref.cpp 
         solid.cpp specpoin.cpp spline3d.cpp surface.cpp triapprox.cpp zrefine.cpp 
-        python_csg.cpp
+        python_csg.cpp splinesurface.cpp
         )
 if(APPLE)
     set_target_properties( csg PROPERTIES SUFFIX ".so")

libsrc/csg/csg.hpp

@@ -24,9 +24,9 @@
 #include "csgeom.hpp"
 #include "csgparser.hpp"
 
-
 #include "triapprox.hpp"
 #include "algprim.hpp"
+#include "splinesurface.hpp"
 #include "brick.hpp"
 #include "spline3d.hpp"
 #include "manifold.hpp"

libsrc/csg/genmesh.cpp

@@ -338,6 +338,20 @@ namespace netgen
 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
 	fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) );
       }
+
+    int bbccnt = 0;
+    for (int k = 0; k < geom.GetNSurf(); k++){
+      auto splinesurf =  dynamic_cast<const SplineSurface*> (geom.GetSurface(k));
+      if (splinesurf)
+	{
+	  for( int i=0;  i< splinesurf->GetNSplines(); i++)
+	    {
+	      string bcname = *splinesurf->GetBCName(i);
+	      if(bcname != "default"){
+		mesh.SetCD2Name(bbccnt,bcname); bbccnt++; }
+	    }
+	}
+    }
     
 
     //!!

libsrc/csg/python_csg.cpp

@@ -233,6 +233,41 @@ DLL_HEADER void ExportCSG()
            }))
     ;
 
+  bp::class_<SplineGeometry<3>,shared_ptr<SplineGeometry<3>>> ("SplineCurve3d")
+    .def ("AddPoint", FunctionPointer
+          ([] (SplineGeometry<3> & self, double x, double y, double z)
+           {
+             self.geompoints.Append (GeomPoint<3> (Point<3> (x,y,z)));
+             return self.geompoints.Size()-1;
+           }))
+    .def ("AddSegment", FunctionPointer
+          ([] (SplineGeometry<3> & self, int i1, int i2)
+           {
+             self.splines.Append (new LineSeg<3> (self.geompoints[i1], self.geompoints[i2]));
+           }))
+    .def ("AddSegment", FunctionPointer
+          ([] (SplineGeometry<3> & self, int i1, int i2, int i3)
+           {
+             self.splines.Append (new SplineSeg3<3> (self.geompoints[i1], self.geompoints[i2], self.geompoints[i3]));
+           }))
+    ;
+
+  bp::class_<SplineSurface, shared_ptr<SplineSurface>,boost::noncopyable> ("SplineSurface")
+    .def("AddPoint", FunctionPointer
+	 ([] (SplineSurface & self, double x, double y, double z, bool hpref)
+	  {
+	    self.AppendPoint(Point<3>(x,y,z),hpref);
+	    return self.GetNP()-1;
+	  }),
+	 (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)
+	  {
+	    auto str = new string(bcname);
+	    self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str);
+	  }),
+	 (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default"))
+    ;
   
 #if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100)
   bp::register_ptr_to_python<shared_ptr<SPSolid>>();
@@ -412,9 +447,19 @@ DLL_HEADER void ExportCSG()
           }),
          (bp::arg("self"), bp::arg("surface"), bp::arg("solid"))
          )
-         
+    .def("AddSplineSurface", FunctionPointer
+	 ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
+	  {
+	    auto planes = surf->CreatePlanes();
+	    auto spsol = make_shared<SPSolid>(new Solid(&*surf));
+	    for(auto plane : (*planes)){
+	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(plane)));
+	    }
+	    spsol->AddSurfaces(self);
+	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf);
+	  }),
+	 (bp::arg("self"), bp::arg("SplineSurface")))
 
-    
     .def("CloseSurfaces", FunctionPointer
          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, bp::list aslices )
           {

libsrc/csg/splinesurface.cpp

@@ -0,0 +1,115 @@
+
+#include <csg.hpp>
+
+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));
+  geompoints.Last().hpref = hpref;
+}
+  
+  void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname)
+  {
+      splines.Append(spline);
+      bcnames.Append(bcname);
+  }
+    
+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));
+      }
+  }
+
+  Array<Plane*>* SplineSurface :: CreatePlanes() const
+  {
+    auto planes = new Array<Plane*>();
+    auto sol = new Solid(new Plane(splines[0]->GetPoint(0),-(splines[0]->GetTangent(0))));
+    for(auto spline : splines)
+      {
+	planes->Append(new Plane(spline->GetPoint(0),-spline->GetTangent(0)));
+      }
+    return planes;
+  }
+
+
+  void SplineSurface :: Print(ostream & str) const
+{
+  str << "SplineSurface " << endl;
+}
+
+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.cpp~

@@ -0,0 +1,92 @@
+
+#include <csg.hpp>
+#include <mystdlib.h>
+#include "splinesurface.hpp"
+
+
+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));
+  geompoints.Last().hpref = hpref;
+}
+
+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 :: Print(ostream & str) const
+{
+  str << "SplineSurface " << endl;
+}
+
+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

@@ -0,0 +1,47 @@
+#ifndef FILE_SPLINESURFACE
+#define FILE_SPLINESURFACE
+
+
+namespace netgen
+{
+  class SplineSurface : public OneSurfacePrimitive
+  {
+  protected:
+    Array<GeomPoint<3>> geompoints;
+    Array<SplineSeg<3>*> splines;
+    Array<string*> bcnames;
+    
+  public:
+    SplineSurface();
+    virtual ~SplineSurface();
+    
+    const Array<SplineSeg<3>*> & GetSplines() const { return splines; }
+    int GetNSplines() const { return splines.Size(); }
+    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]; }
+    int GetNP() const { return geompoints.Size(); }
+    const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; }
+    string* GetBCName(int i) const { return bcnames[i]; }
+    
+    DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false);
+    void AppendSegment(SplineSeg<3>* spline, string* bcname);
+    Array<Plane*>* CreatePlanes() const;
+    
+    virtual double CalcFunctionValue (const Point<3> & point) const;
+    virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const;
+    virtual double HesseNorm () const;
+    virtual Point<3> GetSurfacePoint () const;
+    virtual void CalcSpecialPoints(Array<Point<3>> & pts) const;
+
+    virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const;
+    
+    
+    virtual void Print (ostream & str) const;
+    
+    
+  };
+
+}
+
+#endif

libsrc/csg/splinesurface.hpp~

@@ -0,0 +1,46 @@
+#ifndef FILE_SPLINESURFACE
+#define FILE_SPLINESURFACE
+
+
+namespace netgen
+{
+  class SplineSurface : public OneSurfacePrimitive
+  {
+  protected:
+    Array<GeomPoint<3>> geompoints;
+    Array<SplineSeg<3>*> splines;
+    Array<string*> bcnames;
+    
+  public:
+    SplineSurface();
+    virtual ~SplineSurface();
+    
+    const Array<SplineSeg<3>*> & GetSplines() const { return splines; }
+    int GetNSplines() const { return splines.Size(); }
+    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]; }
+    int GetNP() const { return geompoints.Size(); }
+    const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; }
+    
+    DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false);
+    void AppendSegment(SplineSeg<3>* spline, string* bcname);
+    Array<Plane*>* CreatePlanes() const;
+    
+    virtual double CalcFunctionValue (const Point<3> & point) const;
+    virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const;
+    virtual double HesseNorm () const;
+    virtual Point<3> GetSurfacePoint () const;
+    virtual void CalcSpecialPoints(Array<Point<3>> & pts) const;
+
+    virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const;
+    
+    
+    virtual void Print (ostream & str) const;
+    
+    
+  };
+
+}
+
+#endif

libsrc/include/nginterface.h

@@ -125,6 +125,9 @@ extern "C" {
   DLL_HEADER char * Ng_GetBCNumBCName (int bcnr);
   //void Ng_GetBCNumBCName (int bcnr, char * name);
 
+  // Get BCName for bc-number of co dim 2
+  DLL_HEADER char * Ng_GetCD2NumCD2Name (int cd2nr);
+
   // Get normal vector of surface element node
   // DLL_HEADER void Ng_GetNormalVector (int sei, int locpi, double * nv);     
   

libsrc/interface/nginterface.cpp

@@ -542,6 +542,11 @@ char * Ng_GetBCNumBCName (int bcnr)
   return const_cast<char *>(mesh->GetBCName(bcnr).c_str());
 }
 
+char * Ng_GetCD2NumCD2Name (int cd2nr)
+{
+  return const_cast<char *>(mesh->GetCD2Name(cd2nr).c_str());
+}
+
 
 // Inefficient (but maybe safer) version:
 //void Ng_GetBCNumBCName (int bcnr, char * name)

libsrc/interface/nginterface_v2.cpp

@@ -456,6 +456,22 @@ namespace netgen
         }
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  ElementTransformation<1,3> (int elnr,
+			      const double * xi,
+			      double * x,
+			      double * dxdxi) const
+  {
+    Point<3> xg;
+    Vec<3> dx;
+    mesh->GetCurvedElements().CalcSegmentTransformation(xi[0],elnr,xg,dx);
+    if(x)
+      for(int i=0;i<3;i++) x[i] = xg(i);
+
+    if(dxdxi)
+      for(int i=0;i<3;i++) dxdxi[i] = dx(i);
+  }
+
 
   template <> DLL_HEADER void Ngx_Mesh ::
   ElementTransformation<2,2> (int elnr,

libsrc/meshing/meshclass.cpp

@@ -40,6 +40,7 @@ namespace netgen
     geomtype = NO_GEOM;
 
     bcnames.SetSize(0);
+    cd2names.SetSize(0);
 
 #ifdef PARALLEL
     paralleltop = new ParallelMeshTopology (*this);
@@ -72,6 +73,9 @@ namespace netgen
     for (int i = 0; i < bcnames.Size(); i++ )
       delete bcnames[i];
 
+    for (int i = 0; i < cd2names.Size(); i++)
+      delete cd2names[i];
+
 #ifdef PARALLEL
     delete paralleltop;
 #endif
@@ -94,6 +98,11 @@ namespace netgen
       if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] );
       else bcnames[i] = 0;
 
+    cd2names.SetSize(mesh2.cd2names.Size());
+    for (int i=0; i < mesh2.cd2names.Size(); i++)
+      if (mesh2.cd2names[i]) cd2names[i] = new string(*mesh2.cd2names[i]);
+      else cd2names[i] = 0;
+
     return *this;
   }
 
@@ -126,6 +135,8 @@ namespace netgen
 
     for ( int i = 0; i < bcnames.Size(); i++ )
       if ( bcnames[i] ) delete bcnames[i];
+    for (int i= 0; i< cd2names.Size(); i++)
+      if (cd2names[i]) delete cd2names[i];
 
 #ifdef PARALLEL
     delete paralleltop;
@@ -610,6 +621,17 @@ namespace netgen
           outfile << i+1 << "\t" << GetBCName(i) << endl;
         outfile << endl << endl;
       }
+    int cntcd2names = 0;
+    for (int ii = 0; ii<cd2names.Size(); ii++)
+      if(cd2names[ii]) cntcd2names++;
+
+    if(cntcd2names)
+      {
+	outfile << "\n\ncd2names" << endl << cd2names.Size() << endl;
+	for (i=0; i<cd2names.Size(); i++)
+	  outfile << i+1 << "\t" << GetCD2Name(i) << endl;
+	outfile << endl << endl;
+      }
 
     /*
       if ( GetDimension() == 2 )
@@ -1075,10 +1097,36 @@ namespace netgen
                   }
 
               }
-
-
           }
 
+	if ( strcmp (str, "cd2names" ) == 0)
+	  {
+	    infile >> n;
+	    Array<int,0> cd2nrs(n);
+	    SetNCD2Names(n);
+	    for( i=1; i<=n; i++)
+	      {
+		string nextcd2name;
+		infile >> cd2nrs[i-1] >> nextcd2name;
+		cd2names[cd2nrs[i-1]-1] = new string(nextcd2name);
+	      }
+	    if (GetDimension() == 2)
+	      {
+		throw NgException("co dim 2 elements not implemented for dimension 2");
+	      }
+	    else
+	      {
+		for (i = 1; i<= GetNSeg(); i++)
+		  {
+		    Segment & seg = LineSegment(i);
+		    if ( seg.cd2i <= n )
+		      seg.SetBCName (cd2names[seg.cd2i-1]);
+		    else
+		      seg.SetBCName(0);
+		  }
+	      }
+	  }
+
         if (strcmp (str, "singular_points") == 0)
           {
             infile >> n;
@@ -5714,6 +5762,46 @@ namespace netgen
       return defaultstring;
   }
 
+  void Mesh :: SetNCD2Names( int ncd2n )
+  {
+    if (cd2names.Size())
+      for(int i=0; i<cd2names.Size(); i++)
+	if(cd2names[i]) delete cd2names[i];
+    cd2names.SetSize(ncd2n);
+    cd2names = 0;
+  }
+
+  void Mesh :: SetCD2Name ( int cd2nr, const string & abcname )
+  {
+    if (cd2nr >= cd2names.Size())
+      {
+	int oldsize = cd2names.Size();
+	cd2names.SetSize(cd2nr+1);
+	for(int i= oldsize; i<= cd2nr; i++)
+	  cd2names[i] = nullptr;
+      }
+    if (cd2names[cd2nr]) delete cd2names[cd2nr];
+    if (abcname != "default")
+      cd2names[cd2nr] = new string(abcname);
+    else
+      cd2names[cd2nr] = nullptr;
+  }
+
+  const string & Mesh :: GetCD2Name (int cd2nr) const
+  {
+    static string defaultstring  = "default";
+    if (!cd2names.Size())
+      return defaultstring;
+
+    if (cd2nr < 0 || cd2nr >= cd2names.Size())
+      throw NgException ("illegal bc-number");
+
+    if (cd2names[cd2nr])
+      return *cd2names[cd2nr];
+    else
+      return defaultstring;
+  }
+
   void Mesh :: SetUserData(const char * id, Array<int> & data)
   {
     if(userdata_int.Used(id))

libsrc/meshing/meshclass.hpp

@@ -87,6 +87,9 @@ namespace netgen
     /// labels for boundary conditions
     Array<string*> bcnames;
 
+    /// labels for co dim 2 bboundary conditions
+    Array<string*> cd2names;
+
     /// Periodic surface, close surface, etc. identifications
     Identifications * ident;
 
@@ -588,6 +591,11 @@ namespace netgen
 
     const string & GetBCName ( int bcnr ) const;
 
+    DLL_HEADER void SetNCD2Names (int ncd2n);
+    DLL_HEADER void SetCD2Name (int cd2nr, const string & abcname);
+
+    const string & GetCD2Name (int cd2nr ) const;
+
     string * GetBCNamePtr (int bcnr) const
     { return bcnr < bcnames.Size() ? bcnames[bcnr] : nullptr; }
 

libsrc/meshing/meshtype.hpp

@@ -826,7 +826,9 @@ namespace netgen
     unsigned int seginfo:2;
 
     /// surface decoding index
-    int si;          
+    int si;
+    /// co dim 2 deconding index
+    int cd2i;
     /// domain number inner side
     int domin;
     /// domain number outer side