support of 0D-elements in CSG

libsrc/csg/csgeom.hpp

@@ -115,8 +115,20 @@ namespace netgen
     /// all top level objects: solids and surfaces
     Array<TopLevelObject*> toplevelobjects;
 
+  public:
     /// additional points specified by user
-    Array<Point<3> > userpoints;
+    class UserPoint : public Point<3>
+    {
+      int index;
+    public:
+      UserPoint() = default;
+      UserPoint (Point<3> p, int _index) : Point<3>(p), index(_index) { ; }
+      int GetIndex() const { return index; }
+    };
+    
+  private:
+    // Array<Point<3> > userpoints;
+    Array<UserPoint> userpoints;
     Array<double> userpoints_ref_factor;
 
     mutable Array<Point<3> > identpoints;
@@ -213,11 +225,13 @@ namespace netgen
     void RemoveTopLevelObject (Solid * sol, Surface * surf = NULL); 
 
 
-    void AddUserPoint (const Point<3> & p, double ref_factor = 0)
-    { userpoints.Append (p); userpoints_ref_factor.Append (ref_factor); }
+    void AddUserPoint (const Point<3> & p, double ref_factor = 0)      
+    { userpoints.Append (UserPoint(p,1)); userpoints_ref_factor.Append (ref_factor); }
+    void AddUserPoint (const UserPoint up, double ref_factor = 0)
+    { userpoints.Append (up); userpoints_ref_factor.Append (ref_factor); }
     int GetNUserPoints () const
     { return userpoints.Size(); }
-    const Point<3> & GetUserPoint (int nr) const
+    const UserPoint & GetUserPoint (int nr) const
     { return userpoints[nr]; }
     double GetUserPointRefFactor (int nr) const
     { return userpoints_ref_factor[nr]; }

libsrc/csg/genmesh.cpp

@@ -25,11 +25,14 @@ namespace netgen
     const char * savetask = multithread.task;
     multithread.task = "Find points";
 
+    mesh.pointelements.SetSize(0);
     for (int i = 0; i < geom.GetNUserPoints(); i++)
       {
-	mesh.AddPoint(geom.GetUserPoint (i));
+        auto up = geom.GetUserPoint(i);
+	auto pnum = mesh.AddPoint(up);
 	mesh.Points().Last().Singularity (geom.GetUserPointRefFactor(i));
 	mesh.AddLockedPoint (PointIndex (i+1));
+        mesh.pointelements.Append (Element0d(pnum, up.GetIndex()));
       }
 
     SpecialPointCalculation spc;

libsrc/csg/python_csg.cpp

@@ -532,6 +532,12 @@ DLL_HEADER void ExportCSG(py::module &m)
           }),
          py::arg("solid1"), py::arg("solid2")
          )
+
+    .def("AddPoint", [] (CSGeometry & self, Point<3> p, int index) -> CSGeometry&
+         {
+           self.AddUserPoint(CSGeometry::UserPoint(p, index));
+           return self;
+         })
     
     .def("GetTransparent", FunctionPointer
          ([] (CSGeometry & self, int tlonr)

libsrc/include/nginterface_v2_impl.hpp

@@ -207,6 +207,13 @@ const string &  Ngx_Mesh :: GetMaterialCD<2> (int region_nr) const
   return mesh->GetCD2Name(region_nr);
 }
 
+template <> NGX_INLINE DLL_HEADER
+const string &  Ngx_Mesh :: GetMaterialCD<3> (int region_nr) const
+{
+  static string def("default");
+  return def;
+}
+
 
 
 

libsrc/interface/nginterface_v2.cpp

@@ -475,7 +475,19 @@ namespace netgen
       for(int i=0;i<3;i++) dxdxi[i] = dx(i);
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  ElementTransformation<0,3> (int elnr,
+			      const double * xi,
+			      double * x,
+			      double * dxdxi) const
+  {
+    PointIndex pi = mesh->pointelements[elnr].pnum;
+    Point<3> xg = mesh->Point(pi);
+    if (x)
+      for(int i=0;i<3;i++) x[i] = xg(i);
+  }
 
+  
   template <> DLL_HEADER void Ngx_Mesh ::
   ElementTransformation<2,2> (int elnr,
                                  const double * xi, 
@@ -604,6 +616,16 @@ namespace netgen
     mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
 
+  template <> DLL_HEADER void Ngx_Mesh ::
+  MultiElementTransformation<0,3> (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,3> (elnr, xi+i*sxi, x+i*sx, dxdxi+i*sdxdxi);
+  }
+    
   template <> DLL_HEADER void Ngx_Mesh :: 
   MultiElementTransformation<1,2> (int elnr, int npts,
                                    const double * xi, size_t sxi,
@@ -837,6 +859,30 @@ namespace netgen
     */
   }
 
+  template<> DLL_HEADER void Ngx_Mesh :: 
+  MultiElementTransformation<0,3> (int elnr, int npts,
+                                   const __m256d * xi, size_t sxi,
+                                   __m256d * x, size_t sx,
+                                   __m256d * dxdxi, size_t sdxdxi) const
+  {
+    for (int i = 0; i < npts; i++)
+      {
+        double hxi[4][1];
+        double hx[4][3];
+        for (int j = 0; j < 4; j++)
+          for (int k = 0; k < 1; k++)
+            hxi[j][k] = ((double*)&(xi[k]))[j];
+        MultiElementTransformation<0,3> (elnr, 4, &hxi[0][0], 2, &hx[0][0], 3, (double*)nullptr, 0);
+        for (int j = 0; j < 4; j++)
+          for (int k = 0; k < 3; k++)
+            ((double*)&(x[k]))[j] = hx[j][k];
+        xi += sxi;
+        x += sx;
+        dxdxi += sdxdxi;
+      }
+  }
+
+  
 #endif
   
 

libsrc/meshing/meshclass.hpp

@@ -455,7 +455,7 @@ namespace netgen
     /// Refines mesh and projects points to true surface
     // void Refine (int levels, const CSGeometry * geom);
   
-
+    
     bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const
     {
       if(!boundaryedges)

libsrc/meshing/meshtype.hpp

@@ -1001,6 +1001,9 @@ namespace netgen
   public:
     PointIndex pnum;
     int index;
+    Element0d () = default;
+    Element0d (PointIndex _pnum, int _index)
+      : pnum(_pnum), index(_index) { ; } 
   };
 
   ostream & operator<<(ostream  & s, const Element0d & el);