merge master into branch

libsrc/csg/genmesh.cpp

@@ -692,6 +692,8 @@ namespace netgen
 			     mparam.grading);
 
 	    mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
+            for (auto mspnt : mparam.meshsize_points)
+              mesh -> RestrictLocalH (mspnt.pnt, mspnt.h);
 	  }
 
 	spoints.SetSize(0);

libsrc/csg/revolution.cpp

@@ -156,25 +156,42 @@ namespace netgen
   {
     if(spline_coefficient.Size() == 0)
       spline->GetCoeff(spline_coefficient);
+    if(spline_coefficient_shifted.Size() == 0)
+        spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
 
     Point<2> p;
     CalcProj(point,p);
 
-    return spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1)
+    double val = spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1)
       + spline_coefficient(2)*p(0)*p(1) + spline_coefficient(3)*p(0)
       + spline_coefficient(4)*p(1) + spline_coefficient(5);
+
+    Vec<2> pr = p-spline->StartPI();
+
+
+    // cout << "spline_coefficinet = " << spline_coefficient << endl;
+    // cout << "shifted = " << spline_coefficient_shifted << endl;
+    
+    double val2 = spline_coefficient_shifted(0)*pr(0)*pr(0) + spline_coefficient_shifted(1)*pr(1)*pr(1)
+      + spline_coefficient_shifted(2)*pr(0)*pr(1) + spline_coefficient_shifted(3)*pr(0)
+      + spline_coefficient_shifted(4)*pr(1) + spline_coefficient_shifted(5);
+
+    // cout << "val = " << val << " =?= " << val2 << endl;
+    return val2;
   }
 
   void RevolutionFace :: CalcGradient (const Point<3> & point, Vec<3> & grad) const
   {
     if(spline_coefficient.Size() == 0)
       spline->GetCoeff(spline_coefficient);
+    if(spline_coefficient_shifted.Size() == 0)
+        spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
 
     Vec<3> point_minus_p0 = point-p0;
 
     Point<2> p;
     CalcProj0(point_minus_p0,p);
-
+    /*
     const double dFdxbar = 2.*spline_coefficient(0)*p(0) + spline_coefficient(2)*p(1) + spline_coefficient(3);
 
     if(fabs(p(1)) > 1e-10)
@@ -193,6 +210,27 @@ namespace netgen
 	grad(2) = dFdxbar*v_axis(2);
 	//(*testout) << "grad2("<<point<<") = " << grad << endl;
       }
+    */
+    Vec<2> pr = p-spline->StartPI();
+    const double dFdxbar = 2.*spline_coefficient_shifted(0)*pr(0) + spline_coefficient_shifted(2)*pr(1) + spline_coefficient_shifted(3);
+
+    if(fabs(p(1)) > 1e-10)
+      {
+	const double dFdybar = 2.*spline_coefficient_shifted(1)*pr(1) + spline_coefficient_shifted(2)*pr(0) + spline_coefficient_shifted(4);
+
+	grad(0) = dFdxbar*v_axis(0) + dFdybar * ( point_minus_p0(0)-v_axis(0)*p(0) )/p(1);
+	grad(1) = dFdxbar*v_axis(1) + dFdybar * ( point_minus_p0(1)-v_axis(1)*p(0) )/p(1);
+	grad(2) = dFdxbar*v_axis(2) + dFdybar * ( point_minus_p0(2)-v_axis(2)*p(0) )/p(1);
+	//(*testout) << "grad1("<<point<<") = " << grad << endl;
+      }
+    else
+      {
+	grad(0) = dFdxbar*v_axis(0);
+	grad(1) = dFdxbar*v_axis(1);
+	grad(2) = dFdxbar*v_axis(2);
+	//(*testout) << "grad2("<<point<<") = " << grad << endl;
+      }
+    
   }
 
   
@@ -559,8 +597,17 @@ namespace netgen
 
     if (!spline -> InConvexHull(p2d, eps)) return false;
 
+    /*
     double val = spline_coefficient(0)*p2d(0)*p2d(0) + spline_coefficient(1)*p2d(1)*p2d(1) + spline_coefficient(2)*p2d(0)*p2d(1) +
       spline_coefficient(3)*p2d(0) + spline_coefficient(4)*p2d(1) + spline_coefficient(5);
+    */
+    Vec<2> pr = p2d - spline->StartPI();
+    double val = spline_coefficient_shifted(0)*pr(0)*pr(0)
+      + spline_coefficient_shifted(1)*pr(1)*pr(1)
+      + spline_coefficient_shifted(2)*pr(0)*pr(1) 
+      + spline_coefficient_shifted(3)*pr(0)
+      + spline_coefficient_shifted(4)*pr(1)
+      + spline_coefficient_shifted(5);
     
     return (fabs(val) < eps);
       /*

libsrc/csg/revolution.hpp

@@ -18,7 +18,9 @@ namespace netgen
 
     int id;
 
+    // coefficient for implicizt polynomial
     mutable Vector spline_coefficient;
+    mutable Vector spline_coefficient_shifted;
 
 
     Array < Vec<2>* > checklines_vec;

libsrc/gprim/spline.cpp

@@ -190,6 +190,51 @@ namespace netgen
   }
 
 
+  template<int D>
+  void SplineSeg3<D> :: GetCoeff (Vector & u, Point<D> pref) const
+  {
+    DenseMatrix a(6, 6);
+    DenseMatrix ata(6, 6);
+    Vector f(6);
+
+    u.SetSize(6);
+
+    //  ata.SetSymmetric(1);
+
+    double t = 0;
+    for (int i = 0; i < 5; i++, t += 0.25)
+      {
+	Vec<D> p = GetPoint (t)-pref;
+	a(i, 0) = p(0) * p(0);
+	a(i, 1) = p(1) * p(1);
+	a(i, 2) = p(0) * p(1);
+	a(i, 3) = p(0);
+	a(i, 4) = p(1);
+	a(i, 5) = 1;
+      }
+    a(5, 0) = 1;
+
+    CalcAtA (a, ata);
+
+    u = 0;
+    u(5) = 1;
+    a.MultTrans (u, f);
+    ata.Solve (f, u);
+
+    // the sign
+    Point<D> p0 = GetPoint(0);
+    Vec<D> ht = GetTangent(0);
+    Vec<2> tang(ht(0), ht(1));
+
+    double gradx = u(3);
+    double grady = u(4);
+    // double gradx = 2.*u(0)*p0(0) + u(2)*p0(1) + u(3);
+    // double grady = 2.*u(1)*p0(1) + u(2)*p0(0) + u(4);
+    Vec<2> gradn (grady, -gradx);
+  
+    if (tang * gradn < 0) u *= -1;
+  }
+  
 
 
   template<int D>

libsrc/gprim/spline.hpp

@@ -85,6 +85,7 @@ namespace netgen
     void PrintCoeff (ostream & ost) const;
 
     virtual void GetCoeff (Vector & coeffs) const = 0;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const { ; } 
 
     virtual void GetPoints (int n, Array<Point<D> > & points) const;
 
@@ -138,6 +139,7 @@ namespace netgen
     virtual const GeomPoint<D> & EndPI () const { return p2; }
     ///
     virtual void GetCoeff (Vector & coeffs) const;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
     
     virtual string GetType(void) const {return "line";}
 
@@ -186,6 +188,7 @@ namespace netgen
     virtual const GeomPoint<D> & EndPI () const { return p3; }
     ///
     virtual void GetCoeff (Vector & coeffs) const;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
     
     virtual string GetType(void) const {return "spline3";}
 
@@ -394,6 +397,19 @@ namespace netgen
     coeffs[5] = -dx * p1(1) + dy * p1(0);
   }
 
+  template<int D>
+  void LineSeg<D> :: GetCoeff (Vector & coeffs, Point<D> p) const
+  {
+    coeffs.SetSize(6);
+
+    double dx = p2(0) - p1(0);
+    double dy = p2(1) - p1(1);
+
+    coeffs[0] = coeffs[1] = coeffs[2] = 0;
+    coeffs[3] = -dy;
+    coeffs[4] = dx;
+    coeffs[5] = -dx * (p1(1)-p(1)) + dy * (p1(0)-p(0));
+  }
 
 
   template<int D>

libsrc/include/nginterface_v2_impl.hpp

@@ -70,6 +70,10 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
     ret.index = el.si;
   else
     ret.index = el.edgenr;
+  if (mesh->GetDimension() == 2)
+    ret.mat = mesh->GetBCNamePtr(el.si-1);
+  else
+    ret.mat = nullptr;
   ret.points.num = el.GetNP();
   ret.points.ptr = (int*)&(el[0]);
 

libsrc/meshing/meshtype.hpp

@@ -1134,11 +1134,25 @@ namespace netgen
     MeshingParameters ();
     ///
     MeshingParameters (const MeshingParameters & mp2) = default;
+    MeshingParameters (MeshingParameters && mp2) = default;
     ///
     void Print (ostream & ost) const;
     /// 
     // void CopyFrom(const MeshingParameters & other);
 
+    class MeshSizePoint
+    {
+    public:
+      Point<3> pnt;
+      double h;
+      MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
+      MeshSizePoint () = default;
+      MeshSizePoint (const MeshSizePoint &) = default;
+      MeshSizePoint (MeshSizePoint &&) = default;
+      MeshSizePoint & operator= (const MeshSizePoint &) = default;
+      MeshSizePoint & operator= (MeshSizePoint &&) = default;      
+    };
+    Array<MeshSizePoint> meshsize_points;
     
     void (*render_function)(bool) = NULL;
     void Render(bool blocking = false)

libsrc/meshing/python_mesh.cpp

@@ -609,7 +609,13 @@ DLL_HEADER void ExportNetgenMeshing()
     .add_property("maxh", 
                   FunctionPointer ([](const MP & mp ) { return mp.maxh; }),
                   FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = maxh; }))
-                  
+    .def("RestrictH", FunctionPointer
+         ([](MP & mp, double x, double y, double z, double h)
+          {
+            mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint (Point<3> (x,y,z), h));
+          }),
+         (bp::arg("x"), bp::arg("y"), bp::arg("z"), bp::arg("h"))
+         )
     ;
 
   bp::def("SetTestoutFile", FunctionPointer ([] (const string & filename)

python/csg.py

@@ -14,7 +14,11 @@ def VS (obj):
 
 
 def csg_meshing_func (geom, **args):
+    if "mp" in args:
+        return GenerateMesh (geom, args["mp"])
+    else:
         return GenerateMesh (geom, MeshingParameters (**args))
+#     return GenerateMesh (geom, MeshingParameters (**args))
 
 CSGeometry.GenerateMesh = csg_meshing_func