Merge branch 'codim2integrals' into 'master'

Codim2integrals



See merge request !16

external_dependencies/pybind11

@@ -1 +1 @@
-Subproject commit 2b92a49115657712c7dd924248df2286e6143b8d
+Subproject commit 030d10e826b87f8cdf0816aa36b9a515fb7d064d

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/algprim.hpp

@@ -53,6 +53,7 @@ namespace netgen
   /// A Plane (i.e., the plane and everything behind it).
   class Plane : public QuadraticSurface
   {
+  protected:
     /// a point in the plane
     Point<3> p;
     /// outward normal vector 

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/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,26 @@ namespace netgen
 		      layer,
 		      mesh);
 	  }
+	for(int i=0; i<refedges.Size(); i++)
+	  {
+	    auto splinesurface = dynamic_cast<const SplineSurface*>(geometry.GetSurface(refedges[i].surfnr1));
+	    if(splinesurface)
+	      {
+		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);
+	      }
 		
+	      }
+	    
+	  }
 
 	/*
 	  // not available ...
@@ -1157,7 +1176,23 @@ namespace netgen
 			  refedges.Elem(hi).domout = i;
 		      }
 		    else
+		      {
 		      refedges.Elem(hi).tlosurf = i;
+		      for(int kk = 0; kk < geometry.GetNTopLevelObjects(); kk++)
+			{
+			  auto othersolid = geometry.GetTopLevelObject(kk)->GetSolid();
+			  auto othersurf = geometry.GetTopLevelObject(kk)->GetSurface();
+			  if(!othersurf)
+			    {
+			      if(othersolid->IsIn(edgepoints[0])  &&
+				 othersolid->IsIn(edgepoints[edgepoints.Size()-1]))
+				{
+				  refedges.Elem(hi).domin = kk;
+				  refedges.Elem(hi).domout = kk;
+				}
+			    }
+			}
+		      }
 
 		    if(pre_ok[k-1])
 		      edges_priority[hi-1] = 1;

libsrc/csg/genmesh.cpp

@@ -98,6 +98,12 @@ namespace netgen
 	//(*testout) << " to " << mesh.LineSegment(i).si << endl;
       }
 
+    for(int k = 1; k<=mesh.GetNFD(); k++)
+      {
+	*testout << "face: " << k << endl
+		 << "FD: " << mesh.GetFaceDescriptor(k) << endl;
+      }
+
     if (geom.identifications.Size())
       {
 	PrintMessage (3, "Find Identifications");
@@ -339,7 +345,6 @@ namespace netgen
 	fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) );
       }
 
-
     //!!
     
     for (int k = 1; k <= mesh.GetNFD(); k++)

libsrc/csg/python_csg.cpp

@@ -254,6 +254,61 @@ DLL_HEADER void ExportCSG(py::module &m)
            }))
     ;
 
+  py::class_<SplineGeometry<3>,shared_ptr<SplineGeometry<3>>> (m,"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]));
+           }))
+    ;
+
+  py::class_<SplineSurface, shared_ptr<SplineSurface>> (m, "SplineSurface",
+			 "A surface for co dim 2 integrals on the splines")
+    .def("__init__", FunctionPointer  ([](SplineSurface* instance, shared_ptr<SPSolid> base, py::list cuts)
+	     {
+	       auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive());
+	       auto acuts = new Array<OneSurfacePrimitive*>();
+	       for(int i = 0; i<py::len(cuts);i++)
+		 {
+		   py::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!");
+	       new (instance) SplineSurface(primitive,acuts);
+	     }),py::arg("base"), py::arg("cuts")=py::list())
+    .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;
+	  }),
+	 py::arg("x"),py::arg("y"),py::arg("z"),py::arg("hpref")=false)
+    .def("AddSegment", FunctionPointer
+	 ([] (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,maxh);
+	  }),
+	 py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.)
+    ;
   
   py::class_<SPSolid, shared_ptr<SPSolid>> (m, "Solid")
     .def ("__str__", &ToString<SPSolid>)
@@ -427,8 +482,20 @@ DLL_HEADER void ExportCSG(py::module &m)
           }),
          py::arg("surface"), py::arg("solid")
          )
-         
+    .def("AddSplineSurface", FunctionPointer
-
+	 ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
+	  {
+	    auto cuttings = surf->CreateCuttingSurfaces();
+	    auto spsol = make_shared<SPSolid>(new Solid(&*surf));
+	    for(auto cut : (*cuttings)){
+	      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);
+	  }),
+	  py::arg("SplineSurface"))
 
     .def("CloseSurfaces", FunctionPointer
          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, py::list aslices )

libsrc/csg/splinesurface.cpp

@@ -0,0 +1,67 @@
+
+#include <csg.hpp>
+
+namespace netgen
+{
+void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref)
+{
+  auto pp = Point<3>(p);
+  geompoints.Append(GeomPoint<3>(pp,reffac));
+  geompoints.Last().hpref = hpref;
+}
+  
+  void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh)
+  {
+    splines.Append(spline);
+    bcnames.Append(bcname);
+    maxh.Append(amaxh);
+  }
+
+  string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const
+  {
+    
+    double eps = 1e-5;
+    for(int i=0; i<splines.Size(); i++)
+      {
+	auto pp1 = Point<3>(splines[i]->GetPoint(0));
+	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))
+	  {
+	    return bcnames[i];
+	  }
+      }
+    return new string("default");
+  }
+
+  Array<OneSurfacePrimitive*>* SplineSurface :: CreateCuttingSurfaces() const
+  {
+    auto cuttings = new Array<OneSurfacePrimitive*>();
+    for (auto cut : *cuts)
+      cuttings->Append(cut);
+    for(int i = 0; i<splines.Size(); i++)
+      {
+	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 cuttings;
+  }
+  
+  void SplineSurface :: Print(ostream & str) const
+{
+  str << "SplineSurface " << endl;
+}
+
+}

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,79 @@
+#ifndef FILE_SPLINESURFACE
+#define FILE_SPLINESURFACE
+
+
+namespace netgen
+{
+  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(OneSurfacePrimitive* abaseprimitive, Array<OneSurfacePrimitive*>* acuts) :
+      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]; }
+    int GetNP() const { return geompoints.Size(); }
+    const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; }
+    string* GetBCName(int i) const { return bcnames[i]; }
+    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, double amaxh = -1);
+
+    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;
+    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/include/nginterface_v2_impl.hpp

@@ -13,7 +13,10 @@ NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const
 template <>
 NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
 {
-  return (*mesh)[SegmentIndex(nr)].si;
+  if(mesh->GetDimension()==3)
+    return (*mesh)[SegmentIndex(nr)].edgenr;
+  else
+    return (*mesh)[SegmentIndex(nr)].si;
 }
   
 template <>
@@ -63,7 +66,10 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
 
   Ng_Element ret;
   ret.type = NG_ELEMENT_TYPE(el.GetType());
-  ret.index = el.si;
+  if(mesh->GetDimension()==2)
+    ret.index = el.si;
+  else
+    ret.index = el.edgenr;
   if (mesh->GetDimension() == 2)
     ret.mat = mesh->GetBCNamePtr(el.si-1);
   else

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,
@@ -519,6 +535,17 @@ namespace netgen
     if (dxdxi) dxdxi[0] = dx(0);
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  ElementTransformation<0,2> (int elnr,
+			      const double *xi,
+			      double * x,
+			      double * dxdxi) const
+  {
+    PointIndex pnum = mesh->pointelements[elnr].pnum;
+    if (x)
+      for (int i = 0; i< 2; i++) x[i] = (*mesh)[pnum](i);
+  }
+
 
   template <> DLL_HEADER void Ngx_Mesh :: 
   ElementTransformation<0,1> (int elnr, 
@@ -565,6 +592,15 @@ namespace netgen
     mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  MultiElementTransformation<1,3> (int elnr, int npts,
+				   const double * xi, size_t sxi,
+                                   double * x, size_t sx,
+                                   double * dxdxi, size_t sdxdxi) const
+  {
+    mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
+  }
+
   template <> DLL_HEADER void Ngx_Mesh :: 
   MultiElementTransformation<1,2> (int elnr, int npts,
                                    const double * xi, size_t sxi,
@@ -584,6 +620,16 @@ namespace netgen
       ElementTransformation<1,1> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi);
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  MultiElementTransformation<0,2> (int elnr, int npts,
+                                   const double * xi, size_t sxi,
+                                   double * x, size_t sx,
+                                   double * dxdxi, size_t sdxdxi) const
+  {
+    for (int i = 0; i < npts; i++)
+      ElementTransformation<0,2> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi);
+  }
+
 
   template <> DLL_HEADER void Ngx_Mesh :: 
   MultiElementTransformation<0,1> (int elnr, int npts,
@@ -682,6 +728,15 @@ namespace netgen
     */
   }
 
+  template<> DLL_HEADER void Ngx_Mesh ::
+  MultiElementTransformation<0,2> (int elnr, int npts,
+				   const __m256d *xi, size_t sxi,
+				   __m256d * x, size_t sx,
+				   __m256d * dxdxi, size_t sdxdxi) const
+  {
+    cout << "MultiElementtransformation<0,2> simd not implemented" << endl;
+  }
+  
   template<> DLL_HEADER void Ngx_Mesh :: 
   MultiElementTransformation<0,1> (int elnr, int npts,
                                    const __m256d * xi, size_t sxi,
@@ -691,6 +746,30 @@ namespace netgen
     cout << "multi-eltrafo simd called, 0,1,simd" << endl;
   }
 
+  template<> DLL_HEADER void Ngx_Mesh :: 
+  MultiElementTransformation<1,3> (int elnr, int npts,
+                                   const __m256d * xi, size_t sxi,
+                                   __m256d * x, size_t sx,
+                                   __m256d * dxdxi, size_t sdxdxi) const
+  {
+    double hxi[4][1];
+    double hx[4][3];
+    double hdxdxi[4][3];
+    for (int j = 0; j<4;j++)
+      hxi[j][0] = ((double*)&(xi[0]))[j];
+    MultiElementTransformation<1,3> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 2, &hdxdxi[0][0],4);
+    for(int j=0; j<4; j++)
+      for(int k=0; k<3; k++)
+	((double*)&(x[k]))[j] = hx[j][k];
+    for(int j=0; j< 4; j++)
+      for (int k = 0; k<3; k++)
+	((double*) & (dxdxi[k]))[j] = hdxdxi[j][k];
+
+    xi += sxi;
+    x += sx;
+    dxdxi += sdxdxi;
+  }
+  
   template<> DLL_HEADER void Ngx_Mesh :: 
   MultiElementTransformation<1,2> (int elnr, int npts,
                                    const __m256d * xi, size_t sxi,

libsrc/meshing/curvedelems.cpp

@@ -3686,6 +3686,7 @@ namespace netgen
       }
   }
 
+
   template void CurvedElements :: 
   CalcMultiPointSegmentTransformation<2> (SegmentIndex elnr, int npts,
 					  const double * xi, size_t sxi,

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.edgenr-1]);
+		    else
+		      seg.SetBCName(0);
+		  }
+	      }
+	  }
+
         if (strcmp (str, "singular_points") == 0)
           {
             infile >> n;
@@ -1359,7 +1407,9 @@ namespace netgen
                 if(seg.surfnr2 >= 0)  seg.surfnr2 = seg.surfnr2 + max_surfnr;
                 seg[0] = seg[0] +oldnp;
                 seg[1] = seg[1] +oldnp;
+		*testout << "old edgenr: " << seg.edgenr << endl;
                 seg.edgenr = seg.edgenr + oldne;
+		*testout << "new edgenr: " << seg.edgenr << endl;
                 seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne;
 
                 AddSegment (seg);
@@ -5714,6 +5764,48 @@ 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 )
+  {
+    cd2nr--;
+    (*testout) << "setCD2Name on edge " << cd2nr << " to " << abcname << endl;
+    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())
+      return defaultstring;
+
+    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,12 @@ 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;
+    size_t GetNCD2Names() const { return cd2names.Size(); }
+
     string * GetBCNamePtr (int bcnr) const
     { return bcnr < bcnames.Size() ? bcnames[bcnr] : nullptr; }
 

libsrc/meshing/meshtype.cpp

@@ -2419,10 +2419,6 @@ namespace netgen
   }
 
 
-
-
-
-
   Identifications :: Identifications (Mesh & amesh)
     : mesh(amesh)
   {

libsrc/meshing/meshtype.hpp

@@ -251,6 +251,8 @@ namespace netgen
       singular = 0; 
     }
 
+    void Scale(double factor) { *testout << "before: " << x[0] << endl; x[0] *= factor; x[1] *= factor; x[2] *= factor; *testout << "after: " << x[0] << endl;}
+
     int GetLayer() const { return layer; }
 
     POINTTYPE Type() const { return type; }
@@ -827,6 +829,8 @@ namespace netgen
 
     /// surface decoding index
     int si;
+    /// co dim 2 deconding index
+    int cd2i;
     /// domain number inner side
     int domin;
     /// domain number outer side

libsrc/meshing/python_mesh.cpp

@@ -229,6 +229,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                                      li.append (py::cast(self.surfnr2));
                                      return li;
                                    }))
+    .def_property_readonly("index", FunctionPointer([](const Segment &self) -> size_t
+		  {
+		    return self.edgenr;
+		  }))
     ;
 
 
@@ -548,6 +552,12 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            }
            ))
 
+    .def ("Scale", FunctionPointer([](Mesh & self, double factor)
+				   {
+				     for(auto i = 0; i<self.GetNP();i++)
+				       self.Point(i).Scale(factor);
+				   }))
+                                            
     ;